{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 01 - Introduction To Causality\n",
    "\n",
    "\n",
    "## Why Bother?\n",
    "\n",
    "First and foremost, you might be wondering: what's in it for me? Here's the answer:\n",
    "\n",
    "## Data Science is Not What it Used to Be (or it Finally Is)\n",
    "\n",
    "The Data Scientist role has been labeled [The Sexiest Job of the 21st Century](https://hbr.org/2012/10/data-scientist-the-sexiest-job-of-the-21st-century) by Harvard Business Review. This was no empty statement. For a decade now, Data Scientist has been in the spotlight. AI experts had [salaries that rivaled those of sports superstars](https://www.economist.com/business/2016/04/02/million-dollar-babies). In the search for fame and fortune, hundreds of young professionals jumped into what seemed a frenetic gold rush to get the Data Science title as quickly as possible. Entire new industries sprang around the hype. Miraculous teaching methods promised to make you a Data Scientist without ever showing you a single math formula. Consulting specialists promised millions if your company could unlock the potential of data. AI, or Machine Learning, has been dubbed the new electricity and data, the new oil. \n",
    "\n",
    "During all this time, economists were trying to answer what is the true impact of education on one’s earnings, biostatisticians were trying to understand if saturated fat led to a higher chance of heart attack, and psychologists were trying to understand if words of affirmation actually lead to happier marriages. Meanwhile, we forgot about those who have been doing “old-fashioned” science with data all along. If we were to be completely honest, data science is not a recent field. We're only now becoming aware of it thanks to a massive wave of free publicity from the media.\n",
    "\n",
    "To use a Jim Collins analogy, think about pouring yourself an ice-cold cup of your favorite beer. If you do this the right way, most cups will be beer, but there will be a 1 finger thick layer of foam at the top. This cup is just like Data Science. \n",
    "\n",
    "1. It's the beer. The statistical foundations, the scientific curiosity, the passion for complex problems. All of this was proven very valuable for hundreds of years.\n",
    "2. It's the foam. The fluffy stuff is built on unrealistic expectations that will eventually go away.\n",
    "\n",
    "This foam might come down crashing faster than you think. As The Economist puts it\n",
    "\n",
    "> The same consultants who predict that AI will have a world-altering impact also report that real managers in real companies are finding AI hard to implement, and that enthusiasm for it is cooling. Svetlana Sicular of Gartner, a research firm, says that 2020 could be the year AI falls onto the downslope of her firm’s well-publicised “hype cycle”. Investors are beginning to wake up to bandwagon-jumping: a survey of European AI startups by MMC, a venture-capital fund, found that 40% did not seem to be using any AI at all.\n",
    "\n",
    "What should we, as Data Scientists - or better yet, as “just” Scientists - do during all this craze? If you are smart, you will learn to ignore the foam as a starter. We are in it for the beer. Math and statistics have been helpful forever, and it is unlikely it will stop now. Second, learn what makes your work valuable, not the latest shining tool that no one figured out how to use.\n",
    "\n",
    "Last but not least, remember that there are no shortcuts. Knowledge in Math and Statistics is valuable precisely because they are hard to acquire. If everyone could do it, excess supply would drive its price down. So **toughen up**! Learn them as well as you can. And heck, why not? Have fun along the way as we embark on this quest only **for the Brave and True**.\n",
    "\n",
    "![img](./data/img/intro/tougher-up-cupcake1.jpg)\n",
    "\n",
    "## Answering a Different Kind of Question\n",
    "\n",
    "Machine Learning is currently very good at answering the type of question of the prediction kind. As Ajay Agrawal, Joshua Gans, and Avi Goldfarb put it in the book Prediction Machines, \"the new wave of artificial intelligence does not actually bring us intelligence but instead a critical component of intelligence - prediction\". We can do all sorts of beautiful things with machine learning. The only requirement is that we frame our problems as prediction ones. Want to translate from English to Portuguese? Then build an ML model that predicts Portuguese sentences when given English sentences. Want to recognize faces? Then create an ML model that predicts the presence of a face in a subsection of a picture. Want to build a self-driving car? Then make one ML model to predict the direction of the wheel and the pressure on the brakes and accelerator when presented with images and sensors from the surroundings of a car.\n",
    "\n",
    "However, ML is not a panacea. It can perform wonders under rigid boundaries and still fail miserably if its data deviates a little from what the model is accustomed to. To give another example from Prediction Machines, \"in many industries, low prices are associated with low sales. For example, in the hotel industry, prices are low outside the tourist season, and prices are high when demand is highest and hotels are full. Given that data, a naive prediction might suggest that increasing the price would lead to more rooms sold.\" \n",
    "\n",
    "ML is notoriously bad at this inverse causality type of problem. They require us to answer \"what if\" questions, which economists call counterfactuals. What would happen if I used another price instead of this price I'm currently asking for my merchandise? What would happen if I do a low sugar one instead of this low-fat diet I'm in? If you work in a bank, giving credit, you will have to figure out how changing the customer line changes your revenue. Or, if you work in the local government, you might be asked to figure out how to make the schooling system better. Should you give tablets to every kid because the era of digital knowledge tells you to? Or should you build an old-fashioned library? \n",
    "\n",
    "At the heart of these questions, there is a causal inquiry we wish to know the answer to. Causal questions permeate everyday problems, like figuring out how to make sales go up. Still, they also play an essential role in dilemmas that are very personal and dear to us: do I have to go to an expensive school to be successful in life (does education cause earnings)? Does immigration lower my chances of getting a job (does immigration causes unemployment to go up)? Does money transfer to the poor lower the crime rate? It doesn't matter the field you are in. It is very likely you had or will have to answer some type of causal question. Unfortunately for ML, we can't rely on correlation-type predictions to tackle them.\n",
    "\n",
    "Answering this kind of question is more challenging than most people appreciate. Your parents have probably repeated to you that \"association is not causation\", \"association is not causation\". But actually, explaining why that is the case is a bit more involved. This is what this introduction to causal inference is all about. As for the rest of this book, it will be dedicated to **figuring out how to make association be causation**.\n",
    "\n",
    "## When Association IS Causation\n",
    "\n",
    "Intuitively, we kind of know why the association is not causation. If someone tells you that schools that give tablets to their students perform better than those that don’t, you can quickly point out that it is probably the case that those schools with the tablets are wealthier. As such, they would do better than average even without the tablets. Because of this, we can’t conclude that giving tablets to kids during classes will cause an increase in their academic performance. We can only say that tablets in school are associated with high academic performance, as measured by ENEM (sort of the SAT in Brazil, which stands for National High School Exam):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-03-15T11:20:11.358578Z",
     "start_time": "2023-03-15T11:20:09.259840Z"
    },
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.special import expit\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import style\n",
    "\n",
    "style.use(\"fivethirtyeight\")\n",
    "\n",
    "np.random.seed(123)\n",
    "n = 100\n",
    "tuition = np.random.normal(1000, 300, n).round()\n",
    "tablet = np.random.binomial(1, expit((tuition - tuition.mean()) / tuition.std())).astype(bool)\n",
    "enem_score = np.random.normal(200 - 50 * tablet + 0.7 * tuition, 200)\n",
    "enem_score = (enem_score - enem_score.min()) / enem_score.max()\n",
    "enem_score *= 1000\n",
    "\n",
    "data = pd.DataFrame(dict(enem_score=enem_score, Tuition=tuition, Tablet=tablet))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-03-15T11:20:11.468751Z",
     "start_time": "2023-03-15T11:20:11.360175Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAHwCAYAAAAVRjAeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg6klEQVR4nO3deZhkdX3v8fcHRhQCTNSIIYCAsV3w+sSwReMCEVyiUcbEBYyKxiVRLq6JF7wmMTfgmqjJE0lciBI06LgNkEVWETUKiiGJgDgoCAgCCiLoyPq9f5wzWFNU93RBd1X/pt+v56mn65zzO6e+dbq6PvU759enUlVIktSCzaZdgCRJ82VoSZKaYWhJkpphaEmSmmFoSZKaYWhJkpphaEmLLMkZST447ToGJXlRkls30mbfJJVkx0WqoZI8fzG2PeKxdukf77GTeDwtHkNrmUvy4f6Pefh244g27x6x/gZvPEkumWV7/zKizaEjtveeftmpi/F8NxVz7Oc7blOq6/ljPPb2wCcX6HGfmuSzSX6YZF2SC5P8Q5IHL8T2tXQYWgL4At0byODtgUNt1gGHzPNN4O0jtjf8ifpS4GWDM5LcC3gB8N0x6190SbaYdg1D9uLn+3b3ft7vseE+X9Kq6vtV9bO7u50kfwacCHwH+F3gocAfADcDR9zd7WtpMbQEcHP/BjJ4u3qozX8A5wDvnMf2bhyxvR8NtfkY8MAkvzEw71nAdcDnN/YASV6a5IIkP+s/XZ85eBgryR79J+8fJ7kxydmDj5Xk4CTnJ7kpyeVJjkiyYmD5GUmOTvKXSa4EvtfPf1CSTyX5UZLrkpyc5BHz2CebJXlbkh/0NX0wyZb9Nl/cb2+roef450kuTpLhjVXVNev3LXBNP/vagXmH9vvnp0ku63sdK0fsx/2TnNfvx7OT7D7cZqj9rM8/yb7Asf399T2+D8+xreFeeiV5ZZJjk9zQ1/2GjdSzB/AXwP+tqldW1eer6rtV9aWqehXwh3Ose+Rc+yjJtkk+lOT7/evksiTvGlj+2CRf6mu9Icl/JXnyXPXq7jO0NI7XAk9P8lsLsK0b6IJrsLf1cuCDwJyHl/o3qn8A3go8BNgX+KeB5Q8HzqQLwCcAvw68m/71nuRpwD/SvcE+Ang9cAjw50MP9RzgfsB+wBOS3B/4InA18DjgUcCFwBlJ7reR5/ss4L79er8PPIOuR0q/Hwp49sBz2Ax4MfDBumvXWltHtz93A15Et4/+dqjNZsA7gFcCe/fP61+Hw3Ogpo09//8A/nfffH1v79Vj1v3ndL+7R9J9QHr7Rl5vLwB+Crxr1MKqum6OdTe2j46g68UeAMwAzwUuAEiyOXACcFbfZnfgzX0tWkxV5W0Z34APA7cCNw7dThxqc2p//zjgP4HN+ukCnj/Q9hLgphHbe8NQmzfRvVHeCGxDd0jnZuD+g483S83PBK4Htp1l+bHAf62vccTyLwCrh+a9mu5NbIt++gzgW4PboHtT+srQegG+DbxmjnrP6J/z5gPzXt7vp1/op/8W+OLA8icDtwDbz+N3uGP/e9h3I/vspoHf24v6dfYbaHPv/vfx0n56377NjvN9/nSHgWuer73h104BfzvU5pvAW+fYxr8B/z2Px9ql3/5jx9hHxwMfnqXtvTe2z70tzu2OwyFa1s4CDh6aN9snxsPo3kheRNdbGeW9wFFD83443Kiqzk6yFjiILrROrKqrRhwNG3YK3fmLi5OcApwOfLqqftAv3wP4bFXdPsv6Dwc+PjTv88C9gF+l/zQNnDO0jb2APTIwSKW3Jd0n8bmcXVW3DUx/Cdiif7z/Bt4HfCPJblV1Pl0P9F+r6sqNbHekJL8LvAZ4ELAtXa9qC+CXgSsGmn55/Z2qui7JBXQ9j1HuzvOfr3OHpr9H90FmNmEjPfNZV9z4PjoK+FSSPYHTgM8CJ1XV7f2++iBwUpLT6V4/n6mqC+9KLZo/Dw8KYF1VXTR0u2JUw6r6Lt2htiOSbD3L9q4dsb3ZDtN8AHgF8ELg/fMptqpuBPak+2T8LeCPgIv6w4Z3NNvYZoamM2L+T4babEb35vXIodtD6Hoh49ggmavqPLpDby9Nsh3d4cN57Y87bbg7d/cJusNsz6Q7dPVH/eKNDSiZ6xPDQj7/2dw8NF3M/T51IfCgcQfKzGcfVdVJwAOAI+k+0HwEOL0/NEhVvYzuA9IpwD50HzpmPYemhWFo6a54K91r5/8swLY+Qvcp/Ua6P/55qarbqurMqvozujeOK4Hn9YvPAfbvzwuNch7dm8ygx9MdHvzOHA/7Nbpe2vdGhPI1c6wHsNf6N7veo+neoL89MO99dOH9cuD7dJ/s74rHAj+oqjdV1VlV9S26Q4ijPGr9nSS/SNfjvWCWtvN5/jf329p8lm0stI8AWwGvG7Uwyb1nWW9e+6iqrq2q46rqD4Gn0b1udhtY/o2qeldV/TZwNN3vTovI0BLAFkl+ecRt5KfuqroB+FO6AQyjbD1iW/edZVs/BnYAHjHH4bwNJDkgyWvTjRB8ALAK2Ak4v2/yDrog/GiSPZP8apJnJ3l0v/ytwO8lOSzJg5M8h66n8NdVNfxJf9DfAZsDa5I8Lt0/rD62H4X2mxsp+77Ae5M8rB8I8pfAB6pqsDe3/n+W/hQ4er77Y4QLgfsleUmSByZ5Id1gi2EFvCPJ4/sRgP9E17v851m2O5/nf3H/8xlJ7jdHb3xBVNXXgP8HvCXJe/vnsnOSR6f7v8L3zbLqRvdR/7x+N8lDkszQDaC5Ebg03SjKt/fPf+f+tfU4fv4a1GKZ9kk1b9O90Q16qFluvzTQ5tSh9TajG+wwaiDGqG19Y6jNmzZS01wDMR5Pdx7rGuBnwFq6c20ZaLM3cCrdm/ANdOft9h5YfjBdj+JmuvMmRwIrBpafQTdyb/ixdwY+2j/2TXT/U/YRYNc56j2D7vzfO+nO7d3QT281ou27gdvoBz/M83d4p4EYdKF4Vf/8/43uvGEBu/TLX0Q3AOdJ/X64CfgqsOfANvZlYCDGfJ8/8J7+sW9nloEMfbtRAzGeP9Tm1Lm2MdDu6cDJdCNGf0YXSkcBD+qX78LQQIx57KM/Bb5BF1TX0523emy/bHvg08Dl/X64gu5Q98pp/01v6rf0vwBJS0CS1cCWVfX0adciLUWOHpSWgP7cy+PoBgU8ccrlSEuWoSUtDf9Jd97rHVV1xpRrkZYsDw9KkprRbE/r+uuvN20laRO3cuXKDUYxO+RdktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUjBXTLkCSAFatWjXtEpqxZs2aaZcwNfa0JEnNMLQkSc0wtCRJzUhVTbuGu+T6669vs3BJTRg+x7aczyNN08qVKzM4bU9LktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktSMiYVWktcmOS/JN5Icl+ReSe6T5JQka/uf9x5of3iSi5JcmOTJk6pTkrR0TSS0kuwAvArYs6r+F7A5cCBwGHBaVc0Ap/XTJNmtX/5w4CnAUUk2n0StkqSla5KHB1cAWyZZAWwFXAEcABzTLz8GWNXfPwD4WFXdVFUXAxcBe0+wVknSErRiEg9SVd9L8lfApcA64OSqOjnJ/avqyr7NlUm261fZAfjKwCYu7+eNtHbt2kWqXJI6vs9MzszMzKzLJhJa/bmqA4BdgR8Bn0jy/LlWGTGvZms81xOUpIXg+8zSMKnDg/sDF1fVNVV1C/Bp4DeBq5JsD9D/vLpvfzmw08D6O9IdTpQkLWOTCq1LgUcl2SpJgP2AC4ATgIP7NgcDx/f3TwAOTHLPJLsCM8DZE6pVkrRETeqc1llJPgl8HbgV+E/g/cDWwOokL6ELtmf37c9Lsho4v29/SFXdNolaJUlLV6pmPVW0pF1//fVtFi6pCatWrdpges2aNVOpY7lbuXLlBmMcvCKGJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGSumXYBmt2rVqmmX0Iw1a9ZMuwRJE2BPS5LUDENLktQMQ0uS1AzPaS1hS+08zfA5tqVWn6RNnz0tSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMyYWWkl+Mcknk3wzyQVJHp3kPklOSbK2/3nvgfaHJ7koyYVJnjypOiVJS9cke1p/A3y2qh4K/BpwAXAYcFpVzQCn9dMk2Q04EHg48BTgqCSbT7BWSdISNJHQSrIt8HjgaICqurmqfgQcABzTNzsGWNXfPwD4WFXdVFUXAxcBe0+iVknS0rViQo/zQOAa4ENJfg04B3g1cP+quhKgqq5Msl3ffgfgKwPrX97PG2nt2rWLUrTm5n7XcuLrfXJmZmZmXTap0FoB7A4cWlVnJfkb+kOBs8iIeTVb47meoBaP+13Lia/3pWFS57QuBy6vqrP66U/ShdhVSbYH6H9ePdB+p4H1dwSumFCtkqQlaiKhVVXfBy5L8pB+1n7A+cAJwMH9vIOB4/v7JwAHJrlnkl2BGeDsSdQqSVq6JnV4EOBQ4KNJtgC+A7yYLjRXJ3kJcCnwbICqOi/JarpguxU4pKpum2Ct0iZp1apV0y6hWe67+VuzZs2ibXtioVVV5wJ7jli03yztjwSOXMyaJElt8YoYkqRmGFqSpGZM8pyWpCXm2F23nHYJ2gS84OJ1E3sse1qSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmjBVaSZ6Y5OgkJ/bTeyZ5wuKUJknShuYdWkkOBf4eWAs8vp+9DjhiEeqSJOlOxulpvQbYv6reBtzez/sm8JCFLkqSpFHGCa1tgMv6+9X/vAdw84JWJEnSLMYJrTOBw4bmvQr43MKVI0nS7FaM0fZQ4MQkLwO2SXIh8GPg6YtSmSRJQ+YVWkk2Ax4GPA54BLAz3aHCs6vq9rnWlSRpocwrtKrq9iTHV9U2wNn9TZKkiRrrnFaSRy1aJZIkbcQ457S+C/x7kuPpDg2uH0FIVf3ZQhcmSdKwcUJrS2BNf3/Hgfl156aSJC28eYdWVb14MQuRJGljxulpkWQGOAjYAfgecFxVrV2MwiRJGjbOtQefDpwDPBS4lu7yTV9L8oxFqk2SpA2M09N6C3BAVd1xBYwk+wJ/B5ywsGVJknRn4wx53xH4wtC8L7LhoAxJkhbNOKF1LvD6oXmv6+dLkrToxjk8+Aq6aw++mu7/tHYCfgJ4TkuSNBHjDHn/ZpKHAY8CfgW4Ajirqm5ZrOIkSRo079BK8kjgh1X1xYF5OyW5T1X912IUJ0nSoHHOaX2E7ksfB20BHLtw5UiSNLtxQusBVfWdwRlV9W1glwWtSJKkWYwTWpcn2X1wRj99xcKWJEnSaOOMHnw3cHySdwDfBn4V+GPgyMUoTJKkYeOMHvxAkh8BL6Eb7n4Z8Pqq+uQi1SZJ0gbGumBuVX0C+MQi1SJpwl5w8bpplyCNZZwL5h7U/58WSR6c5PNJTk/y0MUrT5KknxtnIMYRdFd3B/hr4KvAmcBRC12UJEmjjHN48H5VdVWSewGPBZ4F3AL8YFEqkyRpyDihdU2SBwGPAL5aVTcl2QrI4pQ2GatWrZp2Cc1y383fmjVrpl3CSMfuuuW0S9AmYJLnRscJrb+k+xLI24Dn9vP2A7yEkyRpIuZ9TquqPgxsD+xYVaf0s88CDlzfJsljFrQ6SZIGjDvk/adD01cPNfl3YNu7W5QkSaOMFVrz0PT5LYC1T3rjtEvQJmDm5LdMuwRpkzTOkPf5qAXeniRJd1jo0JIkadEYWpKkZix0aDV/TkuStHSNPRAjybbA1oPzquqK/uc2C1SXJEl3Mu/QSrI/8H5gZzbsURWw+QLXJUnSnYxzePBo4C3ASuAeA7ctFqEuSZLuZJzDg/cCPlRVty1WMZIkzWWcnta7gTckcbCFJGkqxulpfQo4CTg8yQZfR1JVD1zQqiRJGmGc0Pok8AXgE4Df0S1JmrhxQmtX4Ner6vbFKkaSpLmMc07reOAJi1WIJEkbM05P657ACUm+AFw1uKCqXrigVUmSNMI4oXVef5MkaSrmHVpV9Rd398GSbA58DfheVf1OkvsAHwd2AS4BnlNV1/VtDwdeAtwGvKqqTrq7jy9JattYF8xN8sQkRyc5sZ/eM8k457leDVwwMH0YcFpVzQCn9dMk2Q04EHg48BTgqD7wJEnL2LxDK8mhwN8Da4HH97PXAUfMc/0dgacBHxyYfQBwTH//GGDVwPyPVdVNVXUxcBGw93xrlSRtmsbpab0G2L+q3gasH/b+TeAh81z/PcAbBtYFuH9VXQnQ/9yun78DcNlAu8v7eZKkZWycgRjb8PMgqf7nPYCbN7Zikt8Brq6qc5LsO4/HGnWpqBoxD4C1a9fOY5PS9Pga1XJyd1/vMzMzsy4bJ7TOpDvndOTAvFcBn5vHuo8BnpHkqXQX3t02yUeAq5JsX1VXJtkeuLpvfzmw08D6OwJXzLbxuZ6gtBT4GtVyspiv93EODx4KPDPJJcA2SS4Eng28bmMrVtXhVbVjVe1CN8Di9Kp6PnACcHDf7GC6f2Cmn39gknsm2RWYAc4eo1ZJ0iZonCHvVybZi25AxAPoDhWefTcv6/Q2YHWSlwCX0oUgVXVektXA+cCtwCF+JYokaZzDg1RVAWcl+er6eUk2Gye4quoM4Iz+/g+B/WZpdyQbHoqUJC1z4wx53z3Jl5P8BLilv93a/5QkadGN09M6BjgR+APgp4tTjiRJsxsntHYG/m9/iFCSpIkbZ/TgZ4AnLVYhkiRtzDg9rXsBn0nyReD7gwv8ahJJ0iSME1rn9zdJkqZirK8mSfJE4CBgu/6rRfYEtl206iRJGnBXrvL+LeBx/ex5X+VdkqS7a5JXeZck6W4ZJ7Tu8lXeJUlaCOOE1vqrvA+a71XeJUm628YZPXgocGKSl/Hzq7z/GHj6olQmSdKQu3KV973oro6xEFd5lyRp3u7KVd7Pxu+2kiRNwTjntCRJmipDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktSMsa6IsRzMnPyWaZcgSZqFPS1JUjMMLUlSMwwtSVIzPKc1ZO2T3jjtErQJ8NyotDjsaUmSmmFoSZKaYWhJkpphaEmSmmFoSZKaYWhJkpphaEmSmmFoSZKaYWhJkpphaEmSmuFlnKRl7AUXr5t2CdJY7GlJkpphaEmSmmFoSZKa4TktaRlZs2bNtEtoxqpVqzaYdt8tDfa0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzXDI+5CZk98y7RIkSbOwpyVJaoahJUlqhqElSWrGsj+n5aVZ5s/L2kiaNntakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGZMJLSS7JTkc0kuSHJeklf38++T5JQka/uf9x5Y5/AkFyW5MMmTJ1GnJGlpm1RP61bg9VX1MOBRwCFJdgMOA06rqhngtH6aftmBwMOBpwBHJdl8QrVKkpaoiYRWVV1ZVV/v798AXADsABwAHNM3OwZY1d8/APhYVd1UVRcDFwF7T6JWSdLSNfFzWkl2AX4dOAu4f1VdCV2wAdv1zXYALhtY7fJ+niRpGVsxyQdLsjXwKeA1VfXjJLM2HTGvZmu8du3aBahO43K/aznx9T45MzMzsy6bWGgluQddYH20qj7dz74qyfZVdWWS7YGr+/mXAzsNrL4jcMVs257rCWrxuN+1nPh6XxomNXowwNHABVX1roFFJwAH9/cPBo4fmH9gknsm2RWYAc6eRK2SpKVrUj2txwAvAP4nybn9vDcCbwNWJ3kJcCnwbICqOi/JauB8upGHh1TVbROqVZK0RE0ktKrqi4w+TwWw3yzrHAkcuWhFSZKa4xUxJEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzVgx7QI0u1WrVk27hDktpfrWrFkz7RIkTYA9LUlSMwwtSVIzDC1JUjM8p7WEeZ5GkjZkT0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS5LUDENLktQMQ0uS1AxDS2O59tpreeMb38h111037VIkLUOGlsayevVqLrjgAj7+8Y9PuxRJy9CSDq0kT0lyYZKLkhw27XqWu2uvvZbTTz+dquL000+3tyVp4pZsaCXZHHgv8NvAbsBBSXabblXL2+rVq7n99tsBuP322+1tSZq4VNW0axgpyaOBN1fVk/vpwwGq6q0A119//dIsfBN20EEHsW7dujumt9xyS4477rgpVqRNyapVq6ZdQjOW09cWrVy5MoPTS7anBewAXDYwfXk/T1Oyzz77sGJF9xVsK1asYJ999plyRZKWm6X8JZAZMW9k72rt2rWLXIoA9tprL0499VQAkrDXXnu576Up2NT/7mZmZmZdtpRD63Jgp4HpHYErRjWc6wlqYe2///6cdNJJ7Lfffuyxxx7TLkdalpbze95SDq2vAjNJdgW+BxwIPG+6Jek5z3kOl156Kc997nOnXYo2McvpPI3uuiU7EAMgyVOB9wCbA/9YVUeuX+ZADEna9A0PxFjSoTUXQ0uSNn0tjR6UJGkDhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZzX4JpCRp+bGnJUlqhqElSWrGimkXoOlLchvwPwOzVlXVJbO0vbGqtp5IYdIiSXJf4LR+8peB24Br+um9q+rmqRSmjfKclsYKIkNLm5okbwZurKq/Gpi3oqpunV5Vmo2HB3UnSbZOclqSryf5nyQHjGizfZIzk5yb5BtJHtfPf1KSL/frfiKJAacmJPlwkncl+Rzw9iRvTvLHA8u/kWSX/v7zk5zdv/7fl2TzadW93BhaAtiy/+M7N8lngJ8Bz6yq3YHfAv46SYbWeR5wUlU9Evg14NwkvwS8Cdi/X/drwOsm9iyku+/BdK/f18/WIMnDgOcCj+lf/7cBvz+Z8uQ5LQGs6//4AEhyD+AtSR4P3A7sANwf+P7AOl8F/rFvu6aqzk2yD7Ab8KU+47YAvjyZpyAtiE9U1W0babMfsAfw1f51viVw9WIXpo6hpVF+H7gfsEdV3ZLkEuBegw2q6sw+1J4GHJvkncB1wClVddCkC5YWyE8G7t/Khkej1v8NBDimqg6fWFW6g4cHNcpK4Oo+sH4L2Hm4QZKd+zYfAI4Gdge+AjwmyYP6NlslefAE65YW0iV0r2uS7A7s2s8/DXhWku36Zffp/x40Afa0NMpHgROTfA04F/jmiDb7An+S5BbgRuCFVXVNkhcBxyW5Z9/uTcC3Fr1iaeF9CnhhknPpDod/C6Cqzk/yJuDkJJsBtwCHAN+dVqHLiUPeJUnN8PCgJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGlrREJan1//M2YtmLknxx0jVJ02ZoSYssyY0Dt9uTrBuYXvRr1vUXgj1isR9HmgT/uVhaZINf5dJfEuulVXXq9CqS2mVPS5qSJHv3X+PyoyRXJvm7JFsMNXtqku8k+UGSd/ZXYBi1rYcmOSXJtUkuTPKcfv7L6a4l+Ya+Z3fiIj8taVEZWtL03Aa8Fvgl4NF0Vw9/5VCbZwJ70l0D7wDgD4Y3kuQXgFOAfwa2Aw4Cjkry8Kp6P91lud5RVVtX1dMX6blIE2FoSVNSVedU1Veq6taqugR4H7DPULO3V9W1VXUp8B66QBr2O8AlVfWhfltfp7tu3rMWsXxpKjynJU1JfwX8d9H1pLai+3s8Z6jZZQP3vwv8yohN7Qz8RpIfDcxbARy7YMVKS4Q9LWl6/p7uCvozVbUt8Ea672oatNPA/QcAV4zYzmXA56vqFwduW1fVK/rlXhVbmwxDS5qebYAfAzcmeSjwihFt/iTJvZPsBLwa+PiINv8CPDjJC5Lco7/t1X8tPMBVwAMX4wlIk2ZoSdPzx8DzgBuADzA6kI6nO2R4LvCvdF+4uYGqugF4EnAgXU/s+8DbgfXfaXY0sFs/SnHNgj4DacL8Pi1JUjPsaUmSmmFoSZKaYWhJkpphaEmSmmFoSZKaYWhJkpphaEmSmmFoSZKa8f8BniZhENe2Ky8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,8))\n",
    "sns.boxplot(y=\"enem_score\", x=\"Tablet\", data=data).set_title('ENEM score by Tablet in Class')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get beyond simple intuition, let's first establish some notation. This will be our everyday language to speak about causality. Think of it as the common tongue we will use to identify other brave and true causal warriors, and that will compose our cry in the many battles to come.\n",
    "\n",
    "Let's call $T_i$ the treatment intake for unit i. \n",
    "\n",
    "$\n",
    "T_i=\\begin{cases}\n",
    "1 \\ \\text{if unit i received the treatment}\\\\\n",
    "0 \\ \\text{otherwise}\\\\\n",
    "\\end{cases}\n",
    "$\n",
    "\n",
    "The treatment here doesn't need to be a medicine or anything from the medical field. Instead, it is just a term we will use to denote some intervention for which we want to know the effect. In our case, the treatment is giving tablets to students. As a side note, you might sometimes see $D$ instead of $T$ to denote the treatment.\n",
    "\n",
    "Now, let's call $Y_i$ the observed outcome variable for unit i.\n",
    "\n",
    "The outcome is our variable of interest. We want to know if the treatment has any influence in it. In our tablet example, it would be the academic performance.\n",
    "\n",
    "Here is where things get interesting. The **fundamental problem of causal inference** is that we can never observe the same unit with and without treatment. It is as if we have two diverging roads and we can only know what lies ahead of the one we take. As in Robert Frost poem:\n",
    "\n",
    ">Two roads diverged in a yellow wood,  \n",
    "And sorry I could not travel both  \n",
    "And be one traveler, long I stood  \n",
    "And looked down one as far as I could  \n",
    "To where it bent in the undergrowth;  \n",
    "\n",
    "\n",
    "To wrap our heads around this, we will talk a lot in term of **potential outcomes**. They are potential because they didn't actually happen. Instead they denote **what would have happened** in the case some treatment was taken. We sometimes call the potential outcome that happened, factual, and the one that didn't happen, counterfactual.\n",
    "\n",
    "As for the notation, we use an additional subscript:\n",
    "\n",
    "$Y_{0i}$ is the potential outcome for unit i without the treatment. \n",
    "\n",
    "$Y_{1i}$ is the potential outcome for **the same** unit i with the treatment.\n",
    "\n",
    "Sometimes you might see potential outcomes represented as functions $Y_i(t)$, so beware. $Y_{0i}$ could be $Y_i(0)$ and $Y_{1i}$ could be $Y_i(1)$. Here, we will use the subscript notation most of the time.\n",
    "\n",
    "![img](./data/img/intro/potential_outcomes.png)\n",
    "\n",
    "Back to our example, $Y_{1i}$ is the academic performance for student i if he or she is in a classroom with tablets. Whether or not this is the case, it doesn't matter for $Y_{1i}$. It is the same regardless. If student i gets the tablet, we can observe $Y_{1i}$. If not, we can observe $Y_{0i}$. Notice how in this last case, $Y_{1i}$ is still defined, we just can't see it. In this case, it is a counterfactual potential outcome.\n",
    "\n",
    "With potential outcomes, we can define the individual treatment effect:\n",
    "\n",
    " $Y_{1i} - Y_{0i}$\n",
    " \n",
    "Of course, due to the fundamental problem of causal inference, we can never know the individual treatment effect because we only observe one of the potential outcomes. For the time being, let's focus on something easier than estimating the individual treatment effect. Instead, lets focus on the **average treatment effect**, which is defined as follows.\n",
    "\n",
    "$ATE = E[Y_1 - Y_0]$\n",
    "\n",
    "where, `E[...]` is the expected value. Another easier quantity to estimate is the **average treatment effect on the treated**:\n",
    "\n",
    "$ATT = E[Y_1 - Y_0 | T=1]$\n",
    "\n",
    "Now, I know we can't see both potential outcomes, but just for the sake of argument, let's suppose we could. Pretend that the causal inference deity is pleased with the many statistical battles we fought and has rewarded us with godlike powers to see the potential alternative outcomes. With that power, say we collect data on 4 schools. We know if they gave tablets to its students and their score on some annual academic tests. Here, tablets are the treatment, so $T=1$ if the school provides tablets to its kids. $Y$ will be the test score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-03-15T11:20:11.478665Z",
     "start_time": "2023-03-15T11:20:11.470021Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>i</th>\n",
       "      <th>Y0</th>\n",
       "      <th>Y1</th>\n",
       "      <th>T</th>\n",
       "      <th>Y</th>\n",
       "      <th>TE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>500</td>\n",
       "      <td>450</td>\n",
       "      <td>0</td>\n",
       "      <td>500</td>\n",
       "      <td>-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>600</td>\n",
       "      <td>600</td>\n",
       "      <td>0</td>\n",
       "      <td>600</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>800</td>\n",
       "      <td>600</td>\n",
       "      <td>1</td>\n",
       "      <td>600</td>\n",
       "      <td>-200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>700</td>\n",
       "      <td>750</td>\n",
       "      <td>1</td>\n",
       "      <td>750</td>\n",
       "      <td>50</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   i   Y0   Y1  T    Y   TE\n",
       "0  1  500  450  0  500  -50\n",
       "1  2  600  600  0  600    0\n",
       "2  3  800  600  1  600 -200\n",
       "3  4  700  750  1  750   50"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(dict(\n",
    "    i= [1,2,3,4],\n",
    "    Y0=[500,600,800,700],\n",
    "    Y1=[450,600,600,750],\n",
    "    T= [0,0,1,1],\n",
    "    Y= [500,600,600,750],\n",
    "    TE=[-50,0,-200,50],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $ATE$ here would be the mean of the last column, that is, of the treatment effect:\n",
    "\n",
    "$ATE=(-50 + 0 - 200 + 50)/4 = -50$\n",
    "\n",
    "This would mean that tablets reduced the academic performance of students, on average, by 50 points. The $ATT$ here would be the mean of the last column when $T=1$:\n",
    "\n",
    "$ATT=(- 200 + 50)/2 = -75$\n",
    "\n",
    "This is saying that, for the schools that were treated, the tablets reduced the academic performance of students, on average, by 75 points. Of course we can never know this. In reality, the table above would look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-03-15T11:20:11.488711Z",
     "start_time": "2023-03-15T11:20:11.480950Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>i</th>\n",
       "      <th>Y0</th>\n",
       "      <th>Y1</th>\n",
       "      <th>T</th>\n",
       "      <th>Y</th>\n",
       "      <th>TE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>500.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>600.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>600</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>600.0</td>\n",
       "      <td>1</td>\n",
       "      <td>600</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>750.0</td>\n",
       "      <td>1</td>\n",
       "      <td>750</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   i     Y0     Y1  T    Y  TE\n",
       "0  1  500.0    NaN  0  500 NaN\n",
       "1  2  600.0    NaN  0  600 NaN\n",
       "2  3    NaN  600.0  1  600 NaN\n",
       "3  4    NaN  750.0  1  750 NaN"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(dict(\n",
    "    i= [1,2,3,4],\n",
    "    Y0=[500,600,np.nan,np.nan],\n",
    "    Y1=[np.nan,np.nan,600,750],\n",
    "    T= [0,0,1,1],\n",
    "    Y= [500,600,600,750],\n",
    "    TE=[np.nan,np.nan,np.nan,np.nan],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is surely not ideal, you might say, but can't I still take the mean of the treated and compare it to the mean of the untreated? In other words, can't I just do $ATE=(600+750)/2  - (500 + 600)/2  = 125$? Well, no! Notice how different the results are. You've just committed the gravest sin of mistaking association for causation. To understand why let's look into the main enemy of causal inference.\n",
    "\n",
    "## Bias\n",
    "\n",
    "Bias is what makes association different from causation. Fortunately, it can be easily understood with our intuition. Let's recap our tablets in the classroom example. When confronted with the claim that schools that give tablets to their kids achieve higher test scores, we can refute it by saying those schools will probably achieve higher test scores anyway, even without the tablets. That is because they probably have more money than the other schools; hence they can pay better teachers, afford better classrooms, etc. In other words, it is the case that treated schools (with tablets) are not comparable with untreated schools. \n",
    "\n",
    "Using potential outcomes notation, we would say that $Y_0$ of the treated is different from the $Y_0$ of the untreated. Remember that the $Y_0$ of the treated **is counterfactual**. We can't observe it, but we can reason about it. In this particular case, we can even leverage our understanding of how the world works to go even further. We can say that, probably, $Y_0$ of the treated is bigger than $Y_0$ of the untreated schools. That is because schools that can afford to give tablets to their kids can also afford other factors that contribute to better test scores. Let this sink in for a moment. It takes some time to get used to talking about potential outcomes. Reread this paragraph and make sure you understand it.\n",
    "\n",
    "With this in mind, we can show with elementary math why it is the case that association is not causation. Association is measured by $E[Y|T=1] - E[Y|T=0]$. In our example, this is the average test score for the schools with tablets minus the average test score for those without them. On the other hand, causation is measured by $E[Y_1 - Y_0]$.\n",
    "\n",
    "Let's take the association measurement and replace the observed outcomes with the potential outcomes to see how they relate. For the treated, the observed outcome is $Y_1$. For the untreated, the observed outcome is $Y_0$.\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = E[Y_1|T=1] - E[Y_0|T=0]\n",
    "$\n",
    "\n",
    "Now, let's add and subtract $E[Y_0|T=1]$. This is a counterfactual outcome. It tells what would have been the outcome of the treated, had they not received the treatment.\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = E[Y_1|T=1] - E[Y_0|T=0] + E[Y_0|T=1] - E[Y_0|T=1]\n",
    "$\n",
    "\n",
    "Finally, we reorder the terms, merge some expectations, and lo and behold:\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = \\underbrace{E[Y_1 - Y_0|T=1]}_{ATT} + \\underbrace{\\{ E[Y_0|T=1] - E[Y_0|T=0] \\}}_{BIAS}\n",
    "$\n",
    "\n",
    "This simple piece of math encompasses all the problems we will encounter in causal questions. I cannot stress how important it is that you understand every aspect of it. If you're ever forced to tattoo something on your arm, this equation should be a good candidate for it. It's something to hold onto dearly and understand what is telling us, like some sacred text that can be interpreted 100 different ways. In fact, let's take a deeper look. Let's break it down into some of its implications. First, this equation tells why the association is not causation. As we can see, the association is equal to the treatment effect on the treated plus a bias term. **The bias is given by how the treated and control group differ before the treatment, in case neither of them has received the treatment**. We can now say precisely why we are suspicious when someone tells us that tablets in the classroom boost academic performance. We think that, in this example, $E[Y_0|T=0] < E[Y_0|T=1]$, that is, schools that can afford to give tablets to their kids are better than those that can't, **regardless of the tablets treatment**.\n",
    "\n",
    "Why does this happen? We will talk more about that once we enter confounding, but for now, you can think of bias arising because many things we can't control are changing together with the treatment. As a result, the treated and untreated schools don't differ only on the tablets. They also differ on the tuition cost, location, teachers...\n",
    "For us to say that tablets in the classroom increase academic performance, we would need for schools with and without them to be, on average, similar to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-03-15T11:20:11.602399Z",
     "start_time": "2023-03-15T11:20:11.490391Z"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAGECAYAAAD6PtKMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAB5bUlEQVR4nO3dd5hcVf3H8ff33qnbS3ohARIg9C699y4gTREUVPghVlBUrIA0RUVRVEAQUKQX6b0oRXoJJUBISG+bzdZp9/z+mEmym53ZkmyZ3f28nmee3Tnn3jtnzs7OfOdUc84hIiIiIlIsvIEugIiIiIhIWwpQRURERKSoKEAVERERkaKiAFVEREREiooCVBEREREpKgpQRURERKSoKEAVERkgZvaUmV0z0OVYk5ldb2aPdeO4T8zs/P4ok4gMLwpQRSSvXJDi8twa8xzzmzznOzP7Qpv7nxS43r/zHHN2nuv9NpfXZeA0nHVSz6tu3bjMN4HPtbnmNWb2VJ7jdgA6/O37ipl9wcyeMbN6M2sys7fN7DIzG9+Lj5E2s1N763oisnYUoIpIZ54Fxq5x22CNY1qAs8xso25c79I81/vCGsfMBr7SNsHMYsDJwKwelr/PmVlkoMuwhh1YXbfb5tKOoX2dd8o5V++cq+vGcYudc03rUNZuM7NrgWuBZ4CDgU2BbwBjgO/2RxlEpP8oQBWRziSdcwvWuC1a45j/Aq8Al3fjeo15rrd8jWNuATYws8+0STsWqAOe7uoBzOx0M3vXzFrNbGmuxW1Cm/ztzOwhM1thZo1m9lLbxzKzU8xsupklzGyOmV1oZqE2+U+Z2bVmdoGZzQfm5tKnmNkdZrbczOrM7BEz26IbdeKZ2SVmtiRXpmvMLJ675pdy1ytZ4zn+1MxmmpmtebFc0LjAObcAWJxLXtYm7YU1u+XXbCFt28VvZj8DTgP2bNMKe2our10Xv5mVm9mfzWxxrv5fNrMD2uRPzp1/nJndZ2bNZvaxmZ3cWQWZ2THAl4FTnHPnO+f+65yb5Zx7wjn3ReCCNsceYmav5P5+i8zsj2ZW2iZ/MzN7OFevTbnXyskrnw/gA3/rQWuziPQBBagi0hu+DRxuZnv3wrUayAapbVtRvwpcA3QaMJjZdsDVwMXAxsBewN/b5G9GtgWuDtgH2IZsF7WXyz8UuA64EdiCbMvcWcBP13io44CRwL7APmY2GngOWATsDuwEvA88ZWYju3i+xwK1ufM+DxxBtqWZXD042ne3e8CXgGtc/+xV/SvgH8DzrG6B/VeBY68DDiTbKr4N8B/g32a2yRrHXUK2jrcEbiUbEE7tpAwnAx86527Jl7mytdfMtgTuJfs33ho4BTiM7GtipX8CS4FdyP6Nv0P29QDZ1ucM8C262dosIn3EOaebbrrp1uEGXA+kgcY1bvetccxjud//CbwGeLn7DvhCm2M/ARJ5rve9NY45H9gxl1cObAIkgdFtH69AmT8L1AMVBfJvBN5YWcY8+c8Ct66R9k2ywxgiuftPAR+0vQbwM+CFNc4z4CPgW52U96ncc/bbpH01V0+luftXAs+1yT8QSAFju/E3nJD7O+y1Zh2vcdw1wFP5/q758vNdC5iSe6xD1jjmVeC63O+Tc8d8p01+KPe3/lonz2M6cG83nu+NwEtrpB0JBMCk3P164NROrpHuLF833XTrn9uqbisRkTxeJNsK1VZzgWPPA94DTiXbkpbPVcAf10hbuuZBzrmXzGwGcCLZAPU+59zCPD3aa3oU+BiYaWaPAk8AdzrnluTytwMecs4FBc7fjI6tg08DMWBD4N1c2itrXGMHYDtrM4EsJw501jII2YAq0+b+f4BI7vHeBP4MvG1mmzrnppNtWb7fOTe/i+v2t01zP59ZI/0ZYOc10l5f+YtzLm1mC8l+ASnE6KL1PGczsn/ztp7Onb8p2THMvwKuyQ1TeIps4PtqN64tIv1IAaqIdKbFOfdhdw50zs2y7Gz+C83s1gKHLevu9YC/AmcCE8l2fXenDI1mtj2wK7AfcAZwmZnt65x7ZeVhXV1mjfuWJ33NiUEe8Djw9TzXq++y4PkfL/ugzr1jZs8Bp5vZJWSHABzVw2u2Faz5GEB4Ha7XlXzBZXKN+47Oh5y9Tzb47I5Cf99ss75zF5jZzcBBZId5/NDMLnPOabkskSKiMagi0psuJvu+8v1euNZNZFsfG8m2jHaLcy7jnHvGOfcTsi2m84GTctmvAPvlxnHm8w6w5xppe5Dt4v+4k4d9mWwANdc59+Eat8WdnAewg5n5be7vTDaA+6hN2p+BL5Lt/l8APNTFNTuzCBi3Rto2XZyTJDt5qDPv5H7usUb67m3y1tZNwBQzOyFfpplVtynDmn+/PckGp9NXJjjnPnbO/dE5dyzwE7JfhFbqznMVkT6mAFVEOhMxszF5bnn72p1zDcCPKbzsT1mea9UWuNYKYDywRSdd8u2Y2ZFm9u3cTP31yLY0TmR1cHIZ2aD3ZjPb3sw2NLPPmdnKLuiLgWPM7Dwz28jMjiM7vvTXzrk1W/3a+gPZoOZuM9s9N1t9NzO7yMx26aLYtcBVZjYtN0nrAuCvrv3yTbfnfv4YuLa79VHAY8DxZnaAmW2ca/We1MU5M4FNcjPgR5hZdM0DnHMfAbcBfzSzA81sEzP7HbA53VvhoSDn3O1kJ7vdYNnVE3Y2s/XMbE8z+xvZeiH3ONua2RW5xz8I+D1ws3NutpmVmdlVZraPma1vZtuQbUmd3ubhZgJ7m9k4MxuxLuUWkbWnAFVEOrM72RbINW95g8qca4EZBfK+n+daBZeOctn1OBt6UN464HCyLYwfkA1ILyQ3JtY59xbZmf0jc4/7OnAO2ZnbOOceILecEfA22Rn+fwR+3tmDOucWkm35XALcSbZL+maygV9XY0VvJ7tywXNkZ+0/AHxvjeu3kp0AFCJbv+viUuB+smNtnyU7BOG2Ls65Fvgf2SXFFpMdG5zP6cDDZFs83yA71OIw59x761hmnHOnkG1B3jv3GO+SHdO8kFwA7Jx7k+wQiD1zj38j2ed6Ru4yaaA693zezV1nIatb2CH75Wo7soFqV63fItJHzDkt8yYiUuxy43rjzrnDB7osIiJ9TZOkRESKWG585e5kl9Daf4CLIyLSLxSgiogUt9fIDqm4zDn31ACXRUSkX6iLX0RERESKiiZJiYiIiEhRGbRd/PX19Wr6FRERERnkKisrOyxdqBZUERERESkqClBFREREpKgoQO0nM2YUWrdcVlIddU111DXVUddUR11THXWP6qlrqqO1owBVRERERIqKAlQRERERKSqDdhZ/Ic45GhsbCYJgoIvSTiwWo76+vtev63keZWVlmHWYACciIiIyKA25ALWxsZFoNEokEhnoorQTjUaJxWK9ft1kMkljYyPl5eW9fm0RERGRgTDkuviDICi64LQvRSKRomstFhEREVkXQy5AFREREZHBTQHqGpYtW8Zuu+3GbrvtxkYbbcS0adNW3U8mk+2OPfPMM7nnnns6XOPZZ5/l+OOP7/Rx3nzzTR555JFeLbuIiIjIUDDkxqCuq5qaGp577jkALr74YsrKyjj77LN7/XHeeustXn/9dQ444IBev7aIiAxdZkYklcRLp0hHoqQ8f6CLJNLr1ILaDTfccAN77703u+66KyeffDLNzc2r8p566ikOPvhgtttuOx566KEO5zY1NXHWWWdx4IEHsvvuu3P//feTTCa5+OKLufPOO9ltt9248847+/PpiIjIIBVLJSid8RbhX38f/6dfJfq3X1FWt4iQ01wEGVoUoHbD4YcfzpNPPsl//vMfNt54Y2688cZVebNnz+b+++/n1ltv5Tvf+Q6tra3tzv31r3/NHnvswcMPP8x9993HT37yE1KpFD/4wQ84+uijee655zj66KP7+ymJiMggE0knCT30L7jg6zD9VVgwB555EL57ErE5H+N5+kiXoUNd/N0wffp0LrroIurr62lsbGTfffddlXfUUUfheR4bbrghkyZN4oMPPmh37hNPPMGDDz7IlVdeiZmRSCSYM2dOfz8FEREZ5CItjXDH3zpmZNLwx18Q/dnVtERL+r9gIn1AAWo3/N///R8333wzW2yxBTfffPOqMapAhwXy17zvnOPvf/87EydObLcO6ssvv9y3hRYRkSHDzGDGdHAu/wHzZuO3NIECVBki1B/QDY2NjYwZM4ZUKsVtt93WLu+ee+4hCAJmzpzJrFmzmDp1arv8fffdl7/85S+43JvKG2+8AUBZWRkNDQ398wRERGTw8zrfMbBA6CoyKClA7YYf/ehH7Lvvvhx11FEdAtApU6Zw6KGHcuyxx3LFFVd02C3q3HPPJZVKsffee7Pzzjvzy1/+EoA99tiD999/X5OkRESkS845mLIZFBpnOmkKQUlZ/xZKpA+ZK9RdUOTq6+vzFry+vp7Kysr+Lk6XWltb+2SrUyje59xTM2bM6PAFQNpTHXVNddQ11VHXirGOwpk00WcegOuvaJ8RicIFf6Vp7CT6+zO9GOup2KiOulZZWdmhe0BjUEVERAaBlB/C2/0gwhtvCffeCEsWwLRt4IBjaCmv7vfgVKQvKUAVEREZJBKhCMnx6xP62o/w0ikykShpxaUyBClAFRERGUScc6TMg3BUM6NkyNIkKREREREpKgpQRURERKSoKEAVERERkaKiMah9oKamhk033XTV/ZtvvpnRo0fnPXb8+PHMnTu3v4omIiIiUvQUoPaBeDzebjtUyK6DKiIiIiJdUxd/P2hsbOTYY49ljz32YJddduH+++/vcMyCBQs4+OCD2W233dh5553573//C8ATTzzB/vvvzx577MEpp5xCY2NjfxdfREREpF8N+xbUWQ0pLny1gfnNGcaW+Jy/bTmTysPrdM2WlhZ22203ACZNmsQNN9zA3/72N0aOHMnSpUvZb7/9OOSQQzBbvXHCbbfdxr777ss555xDJpOhubmZpUuXcvnll3P33XdTWlrKb3/7W6666iq+//3vr1P5RERERIrZsA5QZzWkOOrhpcxsyKxKe3lxkrsPrF2nIHXNLv5UKsUvf/lLXnzxRTzPY/78+SxatKjduNRtt92Wr3/966RSKQ499FC23HJLHnroId5//30OPPDAVdfZYYcd1rpcIiIiUjzMjBAOgoCMHyIIgoEuUtEY1gHqha82tAtOAWY2ZLjw1Qb+umdNrz3OrbfeytKlS3n66acJh8NsscUWHcak7rrrrjzwwAM88sgjnHHGGZx99tlUVVWx9957c+211/ZaWURERGTgxZMt+PNmwYO3QWsz7HkIbtNtaY6XadtahnmAOr85kzd9QYH0tbVixQpGjBhBOBzmmWee4dNPP+1wzOzZsxk3bhynnHIKTU1NvPHGG5xzzjmce+65fPzxx2ywwQY0Nzczb948pkyZ0qvlExERkf4TT7bg3/R7ePqB1YmvP4+Nm0TJT6+iKVY2cIUrEsM6QB1b4udNH1MgfW0dd9xxHHfccey1115sscUWbLTRRh2Oee655/j9739PKBSirKyMP/3pT4wYMYKrrrqK0047jUQiAcD555+vAFVERGSQMjP8BXPaB6crzZuFPXonoSO+SHqYN6IO6wD1/G3LeXlxsl03//rl2YlS62LNdU1ra2u5//77icViBY896aSTOOmkkzrk77nnnjz55JPrVB4REREpDiEzeOT2wgc8ejeR/Y8hHSvtv0IVoWEdoE4qD3P3gbVc+GoDC5ozjOmlWfwiIiIieTkHyUTh/FQSGObNpwzzABWyQWpvTogSERERKSRjBnsdBi8+lf+Az+xNKlbSr2UqRlqoX0RERKSfBEFAsME03PqbdMwsLYNjv0xK4ZlaUEVERET6U3OslNIfXAHPPJgdj5pohR33hqO+SHNZ1UAXrygoQBURERHpZ03REvwDjyWy+0EYjlSslBTW9YnDhAJUERERkQGQCRwtw3y2fiEKUHvZsmXLOOKIIwBYtGgRvu9TW1uLc44nn3ySSCQywCUUERERKW4KUHtZTU0Nzz33HAAXX3wxZWVlnH322bS2thKJREin04RCqnYRERGRQvotUjKzbwOnk13c6y3gS0AJ8C9gMvAJcJxzri53/A+A04AM8A3n3MP9VdbeduaZZ1JRUcE777zDVlttRVlZ2arAFWDnnXfmlltuYdKkSfzrX//iz3/+M8lkku23355f//rX+H7v7mwlIiIiUsz6ZR0DMxsPfAPY3jm3OeADJwDnAY8756YCj+fuY2ab5vI3Aw4C/mhmfRKl2aJ5RK++kNjF3yJ69YXYonl98TB89NFH3HPPPVx00UUFj3n//fe58847efjhh3nuuefwfZ9bb721T8ojIiKDg5kRTSeJtzYTdsFAF0ekX/RnX3MIiJtZimzL6TzgB8BeufwbgKeA7wNHArc45xLATDP7ENgReL43C2SL5hG//By8NkGp/9F0Ws79FW7UuN58KA4//PAuW0Kffvpp3njjDfbee28AWltbGTFiRK+WQ0REBo9YKkFo5vtw53VQtwR/462IHn0qrZUjSJtmfMvQ1S8BqnNurpn9CpgNtACPOOceMbPRzrn5uWPmm9mo3CnjgRfaXGJOLq1XRe68rl1wCuAtmkfkzutInHF+rz5WScnqXSFCoRBBsPpbcGtrKwDOOU488UR++tOf9upji4jI4BNNJwn9+2a4+++rExfMgf88TOwXf6V54obtPktEhpJ+CVDNrJpsq+j6wHLgNjP7Qmen5EkruDHtjBkzVv0ei8WIRqPdKldk6aK86W7Z4lVB47pIp9OkUikymQywOhAdO3Ysjz76KK2trbz55pvMmjWLRCLBTjvtxKmnnsqXv/xlRo4cSV1dHY2NjUycOLHTx1mxYgWLFuV/LoNN27+l5Kc66prqqGuqo64NdB1tXl3ePjhdKZ2GP11A8tuXMHNZff8XbA0DXU+Dgeqoo6lTp3aa319d/PsBM51ziwHM7E5gF2ChmY3NtZ6OBVZGWXOAtlHZBLJDAvJq+yTr6+uJxWLdKpTVjsqfXjOy29foTCgUIhwOr+raX3nNY445hjvuuIP99tuPbbfdlilTphCNRtloo4348Y9/zIknnkgQBITDYX71q191+UesqKjoMogdDGbMmNHlcx3uVEddUx11TXXUtYGuI8/zCL/weOEDPv2YSt8b8L/jQNfTYKA6Wjv9FaDOBnYysxKyXfz7Ai8DTcApwCW5n/fkjr8X+IeZXQGMA6YCL/V2oZJHfxn/o+ntuvmDUeNIHv3lXrn+D37wg1W/t22Rjcfj3HXXXXnPOfroozn66KN75fFFRGQQ87qYx6wxqDKE9dcY1BfN7HbgVSANvAb8BSgDbjWz08gGsZ/LHf+Omd0KTM8df5ZzLtPr5Ro1jpZzf0Xkzuuw5UtxVbUkj/5yr0+QEhER6YkgCGDjLbNBqMszwm2DTcjEtQORDF39NovfOfdTYM3ZPwmyran5jr8IKLwmU2+Va9S4Xp8QJSIisq4SJeVEP/91uOn37TOiMTjzfBLReP7gVWQI0JZGIiIiRSjlh/D2PpzwptvAvTfBkgWw+faw/2dpKa/GKTiVIUwBqoiISJFKhCIkJ2xI6Izz8dIpMpEYaec6WddGhgMzw8yG9DJjQy5A9TyPZDJJJBIZ6KL0i2QyidfVQHoRERm0nHOkzINwVF36w5wPxFoasDkzYfkyWH8jMhU1tETWfeWhYjPkAtSysjIaGxtpaWkZ6KK0s2LFCioqKnr9up7nUVZW1uvXFRERkeLhA/FPP4SLvgktTavSvS12pPQbP6cpEh+4wvWBIRegmhnl5eUDXYwOFi1aNCTWKhUREZH+F2teAb/4OiTbbyRkb70Ed11P+PgzSOXd52hwUt+wiIiISBHzPA/74K0Owekqj91DpLmhfwvVxxSgioiIiBQxM4OFBTfUzAau6VT/FagfKEAVERERKWJBEMBGmxc+oHokLhLtvwL1AwWoIiIiMqSFgNLmFZQuX0xpcwOhQTZU0zlHMHY9GFNgLstJZ5Io6f2J2ANpyE2SEhEREVmpJNmCd/8t8OC/INEKsRJiR3yeYP+jaR5EM99b4mWU/vQq+Msl8Prz2SXHyirgpP8jvc2uZIbYmqgKUEVERGRIimZSeP+4Cp789+rE1ma49a94DfVEjvsaSc8fuAL2gHOOxpIKomf/jFBzE5ZKEsRLSJZUZDdvGGLUxS8iIiJDUri5EZ66P3/mw3cMypnvCT9CU3k1jTWjaY6XD8ngFBSgioiIyFC1fGnh3beCDDTW9295pNsUoIqIiMjQFC/pPD869LYIHSoUoIqIiMiQlCmrLDzzfdIUMqVDa+b7UKIAVURERIak1ngZ/OAKqKptn1E7Cs65lNZoFy2sMmA0i19ERESGJOcczTWjiV1yA97cT2DuJzBxA4Jxk2iJl+GG6ASjoUABqoiIiAxZQRDQHC/Dpm6BbbQlzrlsYKrgtKgpQBUREZEhb1VgKoOCxqCKiIiISFFRgCoiIiIiRUUBqoiIiIgUFQWoIiIiIlJUFKCKiIiISFFRgCoiIiIiRUUBqoiIiIgUFa2DKiIiIjIIhV1ApLkBa26CeJxUvJyEPzRCu6HxLERERGRYMjNCoRCk0wSeRyaTGegi9YuSZAvenX+DR++EdBrMCO+wJ6HTzqUpVjrQxVtnClBFRERk0MmYTyaVIb58IaEn7sJfNA+31U7YdrvSXFpJMIR3jQq7AO/2a+Hh21cnOgcvPQWNK4h/55e0hGMDVr7eoABVREREBpWMecyqTzHx7WeI//lnq9Ltf0/DvyopufAamqpHDdmtTSPNK+Cxu/Pm2fRX8RvqoWZwB6iaJCUiIiKr+J5HvLWZeGsTYYozwFuaDhNrrKP6rxd0zGyox/3pQiLJ1v4vWD+x5ibIpAsfsHgeZtZ/BeoDakEVERERAEoTzdhzD8HDd0Ayib/T3kQP/zzNZVVF02XueR4vzG1h39nTIcg/3tTee4NQSxOJcLSfS9dPovHO86tqB33rsVpQRUREhJLWJrj42/D3K2HhXKhbDA/eCt87mZKGuoEu3ipmxoKmNJZo6fzAITxZKlVSBtvskj9zzASCytr+LVAfUIAqIiIyzHmeh/fRdGzm+x0zG1fA3TcUTXd/EARsVhulZYMtCh80ahxBfPDPZC8kEYrgvvZD3JRN22eMGgc/+h0t8bKBKVgvUhe/iIjIMOc7B0/eV/iA5x8nfOxppGIDH/g459isJsJji+Icv/uhxJ+9v/0BZvDVH5AoKYMgGJhC9oOmeBnx7/8av74OFs2DmhEE1SNpiZcN+u59UIAqIiIiBnQ2XjMcwVE8k26qvRR7bTyad6rPYIOpW1HzwI2wbDFuo83hC18nOXoCwRAOTldqicRhZBwbNX51UDoEglNQgCoiMqx5nsfyjE8GI4yjws8Miw92aS/tIHrgMfDfR/MfsO9RpOJlRRP8OOcY4SUoH1NN85jDCe+4B1ECvEg42/1dJOXsL0Px+WoMqojIMNVEhLs/TbPv3fPY6O8zOez+BTy3xNGitothxzlHZtxk2GX/jpnjJuH2P4p0EQZBUdJUkcDiMZLxElr98JAM1oYjvQuJiAxDCQvxuzfqueyVZavSXl/cyiH3zOEv+4zhc5MjmBu6s6Clo5ZIjPiXvoO//1HwwK2QaIW9DsFttj3NRdR6urai6RTh5gZoWgHxUtKlFSQiMQW0RUoBqojIMLQ87fGrV5flzTv3P4vYd+Ikaj0FqMNNSySObbg5oa9virmAtB/ODvkY5EFcaaIZu/ZyePHJVWn+RltQ+u2LaCqtVJBahNTFLyIyzPi+z8f1KYICn8n1iYBlCY1DHa6cc6Qwkub32XhkM6ORCLOTYd5o8JmXitBCuE8eKxqksRuvbBecAtgHb8Fl5xBLNPfJ48q6UQuqiMgwEwQBJaHO2yfCXvHM2JahxcxYlIlw/IPzeHXR6u1ID55cyp/2Hk0ViV59vHBzI/znkfyZMz/AX1EHI7rYmUn6nVpQRUSGGeccE8t8KiL5PwKmVkWoifRzoWTYqA/CfPbfc9sFpwAPftLEOc8uJmG93JLa3Nj5eqhLFw76feuHIgWoIiLDULWf5l8Hj2PNhtSysMfNB46lytf4U+kbi1odby/N30p6x4cNLE/1crBYUgZeJ+FOzSiNQS1C6uIXERmGzAXsWOvx1ufX558fNPDWkgS7jYtzxPpljAyltBaq9AkzY25TqmC+A1akAkb3YnSSKikjvPO+8J88a7xOmkKmoqb3Hkx6jQJUEZFhyncBY0NJzt2iBEcphiOTSVAkW67LEOScY0JZ4S58zyg49GRtJbwQoVO+jSUS8PIzqzM2nAbfvYTWWMmgX6VgKFKAKiIyzGUy6s6X/jMqCtuMjPHa4tYOecdPLac6FPT6l6SmaAnRM8/ProO6YjmUlpMpraA1Glf3fpHSGFQRERHpNxVemtsPGcfOY9vPnD96wzIu2XUkEZfuk8dNhCI0VtTSOGFDGqtH0aJF+ouaWlBFRESk3zjnGOEnuf2g0dQljYZUQFXUozrsiLneXWJKBi8FqCIiItKvnHOUkqI0Aqxc0kyNmdKGAlQREZFhIOwc4UQTDiNZUkam0FZiIkVAAaqIiMgQ5gEljcvh7r/DS09CJEb8gGNwex5MU7R0oIsnkpcCVBERkSGspHE5fO9kaGpYnXjzH+C/j1By3m9ojilIleKjWfwiIiJDVNgFcOd17YPTHJv5Ad7H7+J1tsuSyADRq1JERGSICida4IUnCx/wxH34WmpJipC6+EVERIYoBxCOFD4gEgUz7aRU5HyDaHMj3oo6CDK4yhoS8XLSZgNdtD6jAFVERGSISsXL8Pc7Cm6/Nv8BBx1LWsFpUQu5gNjH78IVP4KG5QBYJEbs1G+R+sw+JEKdfAEZxPqti9/MqszsdjN7z8zeNbOdzazGzB41sxm5n9Vtjv+BmX1oZu+b2YH9VU4REZGhIu0cbv+jYfzkjpl7HExmzETtplTkYivq4MJvrgpOAUi2wl8uIfzpR9gQbUXtzxbU3wEPOeeONbMIUAL8EHjcOXeJmZ0HnAd838w2BU4ANgPGAY+Z2UbOOW0YLSIi0gPN8TLiP7kK793X4en7IRqHgz9HZtwkWiLxLs+XgRPyDB67EzL5t391t1xN5JzLSISj/VyyvtcvAaqZVQB7AKcCOOeSQNLMjgT2yh12A/AU8H3gSOAW51wCmGlmHwI7As/3R3lFRESGCucczbFSvO12x996Z/A80phaTgcBL52CmR8UzLd5s/EzaRiCAWp/dfFvACwG/mZmr5nZNWZWCox2zs0HyP0clTt+PPBpm/Pn5NJERERkLQRBQMo8Ug4Fp4NEEArD+hsXzHfjJ5EJhfuxRP3H+uNFambbAy8AuzrnXjSz3wErgLOdc1VtjqtzzlWb2VXA8865m3Lp1wIPOOfuWHlsfX39qoLPmDGjz5+DiIhIb4pGo8TjcTKZDI2NjUUXNJaUlOB5Hi0tLWQyGmE3EMyMzarLCJ/7hbzd/Okf/4H3QmWkUqkBKN26mTp16qrfKysrOwyk7a8xqHOAOc65F3P3byc73nShmY11zs03s7HAojbHT2xz/gRgXqGLt32SxWrGjBmDopwDSXXUNdVR11RHXVMdda0v68jMiLc04H04HR57CkaMhj0OJlFeTcof+NawcCZNtLkBXvsP1C+HrXciM348LdGSDsfqtdS1da2jjAsIn38l/PZHUF+XTYzG4NTv4CZvxOQhOou/XwJU59wCM/vUzDZ2zr0P7AtMz91OAS7J/bwnd8q9wD/M7Aqyk6SmAi/1R1lFRET6UmlTPVz8bVg4FxKt2cQ7/kb0zPNhhz0HNEgNZ9JEX/sP/OFnq9dGvf0avI02p/Scy2jKE6RK30qbR8uGmxK99Ea8huXZdVArqkmUlJNmaM7gh/6dxX82cHNuBv/HwJfIjoG91cxOA2YDnwNwzr1jZreSDWDTwFmawS8iUnxWbpMZBMEAl2RwyDiPT1qMj0+9jIhvTHINjL7rz0Te+C9cfRHRTbclVVEzYOWLNi6H3/+0Q7p98DY8cAuho79MurhGIgwLGZddjYF42UAXpd/0W4DqnHsd2D5P1r4Fjr8IuKgvyyQiImsnkk4RaV4B77wCmQxsth2psso+We4mZSGaM0bYg1JLF91Yze5qJszf3mvkpy+tIJ2L50vDYW484nvsX/FX4s/eD6//F2+vIwYk4Pd9H154ovABD91O5KDjSA+jIEkGjnaSEhGRHommk4SfuAduvqpdeviAY/CO+wot4VivPE7afD5t9bjwpSW8tLCVcaVhfrh9DduPDFPK4JoUYma8vDjNj15Y0i69KeU49rkm3jzkNDZ5/hFoWDFgC6+bGSxdVPiA1mbMqaVc+ke/7SQlIiJDQ3jR3A7BKQCP3IH//puruv3Xhed5vLHcsf0/P+GODxv5tCHNiwtaOPLfc/ntmytoKdL2FTMjFAplWyPbaHQhfvHSkrznBA6unm1ktt8Dttl5wGbMZzIZ2H73wgdssAmZITohR4qPAlQREem2MA7u+0fhA+66nnCiZZ0fpy7j89XHF5DJ05t/+SvLWJ4uro8vM2O5i/DMEvjuC/X86s0mPk1GSJANVJOB8cmKwq2+0xuN1mnbk6kd019F7sA5R7DeFBi3XsdMM/jyd0lEe2/nKc/zCPs+oVBxftmQgaVXhYiIdJuXTsOyTrqBly/DW8edbTzPo67V+Ki+cED30sJWjprgF8141GVBhIPvmcN7dclVaT9/aSlX7TWaY9ePEvMCtqiN8uSc5rzn7zzCJ7blAbR4ft78/tISL6P0J1fBv/4Czz4I6TSstyGcdi6Jcev3Sn2bGSXNDdj0V+HlZ2DUeNjrUBIVxbHMlhQHBagiItJtmUiE8FafgffeyH/AptuQjq79GNQWF+LpeSlKuohTwkXUgJrxQlz60rJ2welKZz21kL0mrM+ESIpf7DSC3W+f3eGYqG98cYtaWryO5/c35xyN8XLCp3ybyHFfhSBNEImRiJX22sSt0oY6OP8rULd4deLdNxA96yew3R6kfIUmoi5+ERHpgXTgYM9DoSTPTO5QCI7+Eilbu1ZAz/N4ZmGKEx6cx/t1SbYYkb8V1jfYflS8aFpP69Me10+vL5h/98eNeJ7HRuXGPw4aR01sdf1Mqgjz2NHrMSrUcZeggZQyj6aScprKqmmJxHstOI1m0nDdr9oHpyv96UKiLQ298jgy+ClAFRGRHmkuq4KLrsVtstXqxMkbwYXX0Fo1Yq2vuywd4tznsoHLH95YxkW7jKIk1HFG+2/2GEWVXzwBXcZBa77BsjmLW9KYGVHSHDLe5+Xj1+PlEybx2omTeeboCWxVEQyb2fGhlkZ49T/5M4MA3nm1VybZyeCndnQREemRwDmaascQOecyQq3N4BxBrITW6Lq1ajZnHLMbsuNOZ61IcdFLS7jr8Inc93EDry5qZWJZiO9sW8OkEoiQxvc8QslWnOeT8kMD1qIa9x3bjYrxyqLWvPmHTi5bPTM/CKjxktS0mWs0nPY4sCBYvUNVPs2N/VcYKWoKUEVEpMeccyTC0fYL869jgBjxjIhvJHOtkS8uaOHgu2dz0KQyPjMmzuEblDGtLMA5R0lzI/a/p+E/j0A0RuTQE8msvzEtkd6bZd5d5Zbmt3uOZs/bZxGsUQXbjoqxUWUI5wZ+fGlRCIdg3CSYNytvtrfJVkUzdEMGltrRRUSkKFSFAr6wcUW7tMDBA580ctWbdUws9QmCgNKmFdj5p8O1l2cna73xIvzyW/jX/Zp4ct2XuOop5xyblMNzx0xkpzHZCWKlYeO7W1dx9yFjqTQFpyu5wGHHfzW7bNWatt4Jl0yoi18ABagiIlIkQi7NT3asZbtR7VcBiPjG3YeNpyaUya3DehMsWdDxAs8/hr9gzoDsxFTdsIQdr/gq94ZfYObuCd7dbgUXvvhbxt19NdFUot/LU7TM4OVn8b53OUzbGjwfqkdinzsd2/coXH3dgGzzKsVHXfwiIn3MeSHqMwYu20pIMDA7BQ0GVZbgrkPGMqsxw3/ntzCuNMROY+PU+Gl8lyHc2gRP31/4Ao/dRegrPyDVj7sxRYIM/OOP8PF71Hz8HjVr5If3PIREzcAtwF9MkvEyopEI7i+XYvsfBYecAM0NuGcegntuxP32VnXxC6AAVUSkz5gZizMR/vT6cv75wQp8D06dVsmXN62kxkvqg7iACkuyRTlsXZWddBUECVhZVQ7oLPhM9//s/nCiGZ5/vPABzz+Of8QXB2wL02KSMo/I8V+D2R/DrX9dnRGNwfm/p7W0ovDJMqz0KEA1s/2BE4BRzrnDzWx7oMI590SflE5EZBBbEkTY/fbZzG9aHTRd8NJS/vH+Ch777ESqTV2/nckX0KXiJUR33g+evC//SfseQbq/u4id67xVPKG/c1tN0VLi37scf9ki+HA6VI/Arb8xrSXlZOj/4RlSnLo9BtXMzgb+BMwA9sgltwAX9kG5REQGN8/nhun17YLTlT6qT/HgrCZ8f2C3tRyMUnhw7GlQXtkhz222LZkJG/R7y3Q6EoPNdyh8wE57q/V0DS2ROI1jJtG8x6E0bb4jTSUVCk6lnZ5MkvoWsJ9z7hJg5dfT94CNe7tQIiKD3YqMxz8/KLwrzt/fW0FjRh/Ia6O5vBou/TsceTKMmQCTpsKZ5+O+eSEt0ZJ+L08yHIUvfxcieXa+2m53MiPG9nuZBosgCDTURfLqSRd/OfBp7veVr6YwoPUzRETWYGT3WC8k6puWUVlLgXM0llYSOvrLhA85Aed5JHtxr/hOH9sLsSLj4YAyLyDs0jjnaKkdQ/zX/4Dbr4M3XoCyCjj88wTb7kpLJNbldUWkvZ4EqM8A5wEXtUn7BvBkr5ZIRGQIqPQznLFFFV9/amHe/K9vWUXcC4bVLkK9LY2RjpVm7/RxRZoZS4MIv3uljr9NryeRcRy5QRk/22kE48JpMi6gsaKW8Je+Q6i1JRs0x8u0ZJLIWurJF/izgc+a2SdAuZm9D3wO+E5fFExEZDALgoDDJ5eyw+iOrWf7rVfCjqOiwzp48TyPcChEKDQ4FpNZFoTZ765P+d3rdaxIBiQyjltnNPCZf81ifmr1c0iZT0u8jNZoybD++4qsq269M5iZB0wDdge2ACaR7e5/yTmn/0ARkTwqLcFtB4/lf4sSXPNOPb4ZZ25RyVa1ESqG6e5CnucRb6rH3nkVXv0PjB4Pex5CoqKalFecwarneTw1u5WP61Md8hqSAZe/sozLdqoi5Pp/iSuRoapb7wbOucDM7nHOlQMv5W4iItKFKkty4FiffcaOwADfZYbtvuxmRkn9EvjxV6FuyeqMu64n+s0LcFvtTNorvpUNUnj884MVBfPvndnIj3eooXoABhWHXUCkuQFraoBonFRpOclwVBOPpEd8zwglEwR+iLTnF8Xrp0djUM1sJ+fcC31WGhGRISgIAvzc4icD/7Y/cCKpBPzl0vbBKWTXEb3yp8R+fweNZVUDUrbOeEBFpHD0WRb28Kz//7LxZCv+v/8B9/8T0tnW3dDm2xE+66c0lVYWRZAhxc0zo6SpHp59CF75D1RUET38JNLjJtEaHtjJfT0JUGcBD5rZPWS791e98p1zP+ntgomISN8wM5IWoi7lUZcIKAt7VEccZZbu06Am1NIEb76YPzPIwAdv4W2/Z9GN3fRchrO2qOb2GfmXDTtryyqq/Uy/TngLAf4jt8M9f2+Xbm+/Ahd9k9hPrhqQJbdkcCmpXwo//BI01K9O/N/ThI74PNEjvkginGfptH7Skw6JOHA32cB0AjAxd5vQ+8USERleMl6IFUGYZsJ4Xt/2FZePX5/vvbCcaTfNZKdbZ7H5zTP53EMLWZyJYNaHa7N2ttsSQHNj3z32OnDOMbXS48ubddwc4DNjYhw3pbzfg+po8wq496b8mXNm4i9Z0K/lkcEnmknB9b9pH5yudO/NhBvq+r9QbXS7BdU596W+LIiIyHDkzGNBKsRvXqvjsU+bqI35nLtdDbuMjlLWB8tMp70QV7zdyPXT238o/Xd+C0fcN5eHjhjXZxO4gmgcf/R4WDg3/wHTti661tOVykhxwY7VnL5pFddMX05zyvHFaRVsVh2mcgC2rLXWFki0Fj5gzky8iVOKtj5l4IUSLfDqc4UPeOlp/ENPGrBd0Ho0ZdLMpgInAuOBucA/nXMz+qJgIiJDnZnxcYvPrrd9Qks627X+cX2K4x6YxynTKvjlTjWU0nHm+LpYlvK4bnqeFhNg+rIE81sCKvqoZzhRUk7JV74PF36jY+Zn9iZdUd03D9xLSkmxaRlcuUsVDggyfT/hzfezE1bWDDRdNIaFQpAusHLAqHEagyqdcy57KyTdu+89PdXtfiQzOxx4BdgEWEZ2i9OXzeyIPiqbiMiQ1uRCnP3UwlXBaVs3vLuChZ00kK2txlRAqpNGtY/rU33WzR8EAYkNpsFF1+I23hI8D2pGwpe+Q3D69wZ8UkZ3pdNpMum+Ha8bT7ZQNv8T4g/dSsnzj1HWuJxQmyESyZIy2OOQ/CdX1RKMmagAVTqVjsZxm21X+IAd9xqw1lPoWQvqL4EjnXOrdo4ys72APwD39m6xRESGvsaM8Z/5LQXz7/+kkW9tVtKrHxJlYQ/fIFMgdlmvPIzrw/U8U36Y9MQpRM+9HC+TwmEkS8rJqCt6ldJEM/zmh/Du66sTPZ/Yty+idfPtSXshUuYTOeEMrG4JvPbf1ceNGA3n/57W0oo+311LBrdkOEr4tHPhvFMhuca34V0PIF09YkDKtVJPAtQJwLNrpD2HJkmJiPSJUB+0ZFaHAk7YqJyb3+84I32DyjATSvt+MU/nHK3hKKycIaxAapUQDnvglvbBKWQnmP3mh8SuvJ3G8hoAmqIlxM76KaGmelg4DyqqCKpH0lrS/5O2ZPBxztFSO4b4Ff+Ae26EN16A8io48mQy07YZ8B6NngSorwPfBS5tk/adXLqIiPRQhR+w/3qlPDq7KW/+IZNLyWR6dxxY2KX5ybblLE8G3D9z9eNuUh3h7sPGU+WlOh2WJn0r0tIED92ePzMI4H/P4O9/zKpW9dZwFKpGYdWjV3fpKziVbsoAjeU1hD//DULHno7zQySj8aL4gtOTAPVM4D4z+ybZdVAnAk2AxqCKiKyFuGX4zR6j2PW2WdQn2n8gfGebGkZG+yZSbJjzEdfsNY1lu45iflOampjPiKhR6aU0bnGAWZCB1ubCByxZkHeMsP5usi5SZqRipdk7RRCcQs+WmXrPzKYBOwHjgHnAi865gZ3mJSIySDnnmBhJ8fIJk7nx3XoenNXEqJIQ39mmmk0qfUp6eQZ/28eNk2J8GMZXASt3uVKMM+CCcAR/vQ1h9kf5D9hqpwGduCLSX7odoJrZ1sBS59xzbdImmlmNc+6NviiciMhQ55xjpJfgu1uUcsamZUR8CLsM+u4/PCVipZR86bvw8//rmDlmIsGkqWotlWGhJ6PhbwLCa6RFgBt7rzgiIsOTCzKUWJpQ0LfLF0lxC4KAxKSpcP6VMGZiNtHzYdf94ad/pCVeNrAFFOknPRmDup5z7uO2Cc65j8xscu8WSUREZPhK+WEyG29N9Bd/xku0gB8iFS8j6fkahyHDRk8C1Dlmtq1z7tWVCWa2LdmxqCIiItJLgiCgJVoC0T7a1kukyPUkQP0NcI+ZXQZ8BGwInANc1BcFExERGUzCOLxMmkw4QjpQS6fIuujJLP6/mtly4DSyS0x9CnzXOVdgwTYREZGhL5pOEl48D+67GZYtJrzVTrDHwTSXVRGoS15krfSkBRXn3G3AbX1UFhERkUElkkkRfvJeuOkPqxPffR3uvYmSX15LU82YLie9mRm+7wOQTvfdNrMig0m3Z/Gb2Ym5dVAxs43M7Gkze8LMNum74omIiBSvSFND++B0peZG3J8vJrLmHudrWOEivLDM+NqzdXzrv8t5u9Gn0a25YI7I8NOTFtQLgV1yv/8a+B/QCPwR2KeXyyUiIjJgwuEwzYRJBIZvjio/02H7R8/zYPorBa9h775OqLWZRDiaN7/eRfnio/N5em7LqrRr3qnnlGkVXLRTLWUkO17TjFCQwcukSUeiZDTWVYaongSoI51zC80sBuwGHAukgCV9UjIREZEBkMCnrmoSZz+8iNcWtzKxPMQPt69l7/FRytbc3aurLnmXf9tIz/O476OmdsHpSje8u4IvbFLJjtXWbnhALJUgtGB2dqzr8jrC2+0Kux/E/JIRLE86Yr5RHQoIOQ0TkMGvJwHqYjObAmwB/M85lzCzEqDjpsAiImtp5Vg8befY+1bWrRRmnsezCzIce//qFRTfXZbk5Efm842tq/nhNhXEyAaAQRDAptsWvtikqQSx/MtELc/4/P6NuoKnXvlGHdfvXYufCzaj6SSh+/8Bd92w+qD334B7byLzo7+yx1MBjamAkzau5OefqaXKEj141iLFpyc7SV0AvAJcC1yeS9sX0DanIrLOGonwblOIK95q5q/vtzIvFSHRs3mcUkCDC/N2o88lbzRx1bsthMdvRKtT3eazLB3i608uzJv3+9frqEu3/9hMVVTB/kd3PDgUgjN+SCJWmvdagYPGVP7WVYAVyYBMm977cMPy9sHpSo0rGHvrlXxz4yipAG54t54TH5pPg8axyiDXk2WmrjezW3O/N+eSXwROWHmMme3qnPtP7xZRRIa6BqJ87cmFPPBJ06q07zwLv9p9FF+YEl/VYiU9t4IoX3p8AY9/2rwq7Tzgst1GcvLUEuKq23aWpxwLmvPXiQPeXJJg3FhvVdd7IhTFO/6r+Ft/Bu68AeqXwrRt4Jgv01o1ssO41ZXKQwGHr1/G1W8tz5t//NRy4p4jk8m1fL/ybMEy+6//l8OP/jY/y91/fkEL81ocG+dpvPU9j0hzA4Yj44dJRmLaWleKUk+XmWpe4/6iNQ55EKhY10KJyPDheR63zWhsF5yudM6zi9h7wmQ2jA1AwYYAz/f5x7sN7YLTlb733GL2nai6XVOoi0FrJaGOHY8t4RjeljsT2XhrLJMmHYmR8jofTuEHGb69TQ3/eH8FK5Ltg9iJ5SEOWq+ETKbNJKnOxro6RzZ8Xu3NJQmmTQ63C5BLki14/3s62xJbtxibshnhL36DxJhJpHy1qEtx6UkXf3doPKqI9Ehdxue3rxUei/fXt5fjh/ThuTaWpT2ufL1w3V77Tr3Gpa6hOmJsMSL/rPuob0yrCedtcQyCgNZwlJZYaZfB6UpjQklePH4Sx29UTtiD0rBx5pZVPHn0etT6qydjZTIZ2G63gtdxm2zNM8vbP+a40lC7ckZTCbwbr4S/XAKL50M6jb33BvzwNKIfvZNdkUCkiPT2K1L9BCLSI4EzFrcUbh2a15RGK+msHeeMZa2FJ5vNa0rjTO0KbZVbihv2H0t5pP3HowE3HDCWar/3hkQ45xgbSnLVbjV88MX1eefzk7lw+wpGeIkOQXC6qhZ2P6jjRSJR5p/0XS59d3VAWxHx2KS6fSAdblwOzzyYvyB/uYRYc0NvPCWRXqNmCREZUCW+Y8/xJTw0q2MXP8AR65dhLtC337UQ9xx7Tijh4QJ1e9QGZVBgjORw5Zxjw5IMLxwzjn9/muDpOS1sVB3mtE2rGB3J4LveX10i5NJUr/yeUOCF3hqOEf/iN/F32APuuRFWLMdtsQOLDjyZk1+DBc3Z4QDlEY8Hj5pIrZ9edS3P8+CDtwsXYPF8vJYmKDChS2QgKEAVkQEVI82Fu4zk0dlN7WYtA4wtDbH3hDhB0HHBcula3NJctPNIHstTtxPLQ+wxPk4QaDmiNbkgoPnTGfzfJpvwlY1LCJtrPx50gLRE4njb7EZk022xTIZMrISMi3LBbhleW9zKemVhthwRyQWnq794OOcg3nnw6TQGVYqMxqCKyIByzjE5luHpYyex9cjsjB3P4LMblvHUMRMZ4ae6uIIU4pxjcjzDs5+bxPajs3Ub8uD4qWXZcY7ewAddxSyTyeAF6aJakzc71jVGS6yUJEaFJdmyPMOXpkTZZ7Rl/6ZrbA7gnIMNp0GBINRtug3pkrL+KL5It/X4K5OZVQDtXsnOuXm5n+W9VC4RGUZCBGxR7vj3YWNoShu+BxV+QMQl0Qo46yZEwGZljnsOGUNTxvCAVN3C3DjHgS6d9JZCy1mtlCyrJPKzP8KfL4Y5M1dnVFRjZ/yIZDiKXhBSTLodoJrZfsBfgEm0byl1gKaBisg6cc5RRoqyle9K+qzsNc45SklRmnunnlG3BEZUD2yhpF8481icDvHw7GbeXTqeXc/8K7uWtDLh0RthvSm4rXeipbQSp7HIUmR60oJ6LdndpG4BOm4eLCIiIsXDPF6vNw66+xMSuUHIf3wLRsR9njrmm0yKZbLDFxScShHqyRjUGPA351yjcy7T9tZXhRMREZG1sywT4rP/nrMqOF1pSUuGzz88j/q01j6V4tWTV+dvgO+ZadE8ERGRYjerMcPyRP7W0TcWJ6jT/EMpYj3p4r8DeBj4gZktaZvhnNugV0slIiIia83MqE903sGZXHPtsV5+/GiiBT+T3QwiUVJGRjtuSA/0JEC9HXgWuA2NQRURkR4wz2dZ2iftIOI5qv1MlzPPZe0559i4OlIwvyLiURXtmy7+SCZF5NMP4frfwsz3scpq4kd+kWD3g2iOxPvkMWXo6UmAuj6wjXNurd9RzMwHXgbmOucOM7Ma4F/AZOAT4DjnXF3u2B8ApwEZ4BvOuYfX9nFFRGTgNBDh+ukNXPHqMuoSARtUhrl4l5HsPiZMCb23dai0VxUK+PKmlVw3vb5D3gU7j6A2lIFe/o7gex6R6W/ApeesTqyvg7//Du+Dt4id/n1aw9HefVAZknry9ekeYJ91fLxvAu+2uX8e8LhzbirweO4+ZrYpcAKwGXAQ8MdccCsiIoNIM2G+89wSfvz8EiqjPjuMjtGQDDj+wXncO6sV8/TW3lfipPnFZ2q4fLdRjIxn63mDyjD/PGgcn9ugBILen+McbW6Aay/Pn/nCE4Qa6nr9MWVo6kkLahS418yeBRa2zXDOfbGrk81sAnAocBHwnVzykcBeud9vAJ4Cvp9Lv8U5lwBmmtmHwI7A8z0or4iIDLClSWN2fZJX948xet77RBbPpWWLacwoHcdZry3jgIkTqPG0GExfKSPJ1zaJ8bkp62WHV5ijys8QBH0zQ8praYIlCwsf8MHbeLuM0/AO6VJPAtR3cre19Vvge0Db3aZGO+fmAzjn5pvZqFz6eOCFNsfNyaWJiLTjeR7Ouex2jlJUzIy5DQnu2GQFYy74OrQ2A1ANjBszgXu+9TsaUo4a9fj2qSCTocoyq7bY6cvY0IXCne95Xlqu/1XpFuuPF4qZHQYc4pz7PzPbCzgnNwZ1uXOuqs1xdc65ajO7CnjeOXdTLv1a4AHn3B0rj62vr19V8BkzZvT5cxCR4hEKhYiNnMCCVIj36lJMrggxqcSRXjKXRCIBgO/7hMNh0uk06bTGOQ4E3/eZNGIklT/4PDSu6JAfbL49y77+Sz6dM2sASid9YXxVBaOuvxzefrljZihE4lf/ZPqS5f1eLik+U6dOXfV7ZWVlh+81PWlBxcz2Jzs2dJRz7nAz2x6ocM490cWpuwJHmNkhZBf8rzCzm4CFZjY213o6FliUO34OMLHN+ROAeYUu3vZJFqsZM2YMinIOJNVR11RH2Va5xZkIh9w7hw/qkqvSx5SEeOSzUyhbOov4mEl8WJ/m3bokm1RHmFIZosZLqeUmpz9fRyWfvJ83OAXw3n6ZmnQTsSJ8Tet/rXvWrCczgzN+CD/+KtS1WZHS8+A7F+OqaphaPXIASjpw9FpaO90OUM3sbLKTnK4Bjs0ltwBXArt0dq5z7gfAD3LX2YtsC+oXzOxy4BTgktzPe3Kn3Av8w8yuAMYBU4GXultWERm6mlyIUx+d3y44BVjQnOawe+fwyFGT2O32T5nTuLrVdFxpiEc/O5EJkd4NUj3PY3naxwEVIYcFvd9S63ke0ZZGLAjIhKMkw5GiCLR9PzvpJgiCguUxM7z6pZ1ex0u0th/4JX0u5PuYCwjMy2512oucczRV1FJy8d+w996A15+HsevBLvvRWlZFWvOdpZt60oL6LWBf59wnZvb9XNp7wMbr8PiXALea2WnAbOBzAM65d8zsVmA6kAbO0paqIgJQlzKem5d/KebZDSneX55mYXP7QHFeU5qj75/Lo0eOp8KSec/tqXoX5db3G/jL2/UkMo5jp5Tzf1tWMSqU6rUJICXJVrzX/gN33wB1S/E22ZLwSWfROmIs6QGa/Z7EZ2na574ZjSxqznDgpFI2qgxTmadenXMwsZN9XOKluJKyPiyttBXOpIiuqIPH74aF82CrnXDb7kJLaSVBL37pcc7RFC/H224PvB32wjnX64GwDH09CVDLgU9zv698JYeBHr3bO+eeIjtbH+fcUmDfAsddRHbGv4jIKs3pzoO/RS1pSsNehy0e369LsjTpqFjHCTme57E0E+az98/l9cWJVelXvLaMG9+r57/HTWKkl+jkCt0TSyXw/v5bePah1YmvPQ9vvEjsp3+ieYNp/T4TOonPvz9NctpjC1Z9CPzq1WVsOSLKvYeNp8o6Pu90eTWhbXaB1/7b8YLHfYVEafnqTxTpM6EgQ/SNF+DKn8DKYPSlp7DyKkouvIam6pG93jIfBIFm68ta68k6qM+QW6e0jW8AT/ZecUREOlcZ8Qh38s41rjTEimT+D8XGAund1UCEF5bCw7Nb2gWnKy1uyfD71+twXo+G9+cVaqhrH5yuFATwl4uJtjSu82P01JK03y44XenNJQkuf20ZQZ5W3dZwFPd/P4bDToRILJtYWQNfPY/07geTVnDaL2ItjfCHn68OTldqWI770wVEkq0DUzCRAnryLno2cJ+ZfQUoN7P3gRXA4X1SMhGRPGpCAf+3ZTW/e73jgt8HTy7llYWt5Nvy2yC3WPnadTU2uDDfem4J1TG/3fjWNf1rxgq+s3UlVeuwi6TnefDOK4UPmPsJXmszREvW/kF6yPd97v2gsWBj53Xv1PPdravzrmnaFIkT/txXiRx6IqRSuEiURGkFGbWu9Qszgw+nQyb/69bee4NQSxMJ7fAkRaTbb6G59Up3AI4HTiI7qekzzrkFfVQ2EZEOQi7NOdtU8cMdaikNZ1cmCXtw2maV/GnvMcxtzN+9fvImFVSH1i4gMjNmNgbc8WEDmYBOW3DDnnW+DmR3RQvsWV5Vi33l+3jJBGUz36Vs+WJiqb5v/TKzDmN722pJOzKdtIamMJpKKmiqrKU5XqbgtB+ZGSSaOz+mQPAqMlB61A/lsgNUXjSz/61MMzPPOad3GhHpN+UkOWfLMr68aQVNKUdJyKgOB4SCVr6zWQmV0RB/eKOO5rQjHjLO3KKKb21dRbRnQ+ZX8X2fG95dDsAjsxv5xc4jue/j/F3sp29WRVUoYF3eFYMggGnbgFn7LtnqEXjf/AXBdb+Gv166unybb0/p2T+nKVa69g/ahUwmw0GTyvjNa/m3qtx6ZIy4p/76YhQEAUzdovABo8YRxPvutSOyNrrdgmpm25rZ82bWBKRyt3Tup4hIv/KCNCO8JJOiKUb6SUK5JZ5WzP2Y87Yu583PT+atz0/mrZMmc/625VSsZXC60sphA3Mb07SmHQdO6viBvnF1hFOmVeB6YY/zREk5fLX9sH877isE11wOsz9qn/72y9g1lxHN9N3bsXOOTapCbF7bsRvYgN/sMYpyT61wxSpVVgl7Htoxwwy+eh6tJVrrS4pLT1pQbwDuA74MdN5XICKyjlau/+kFAYHnkYiXdXtGcDZ4BSK5hHXs48lkMpy8SQV/m14PwLefWchv9xzNsVMruH3GChIZxxc3qWTvCXFqvGSHeShrI+WHsM/sQ2TjLeHhO2DZImziBrg5M/Of8PIzhJsbSJTXrPuDF1Dlpbjv8PFc+soyrp9eT2vGseWIKL/dYxSbVVAU67NKfolwFO8LZ+NP2xruuh7qluCmbIadfDaJ0RM1216KTk8C1EnAj5zegUSkj8VSCULvvQ7/+CPMm4U/ZgIlx3+NzBY70jIAEzmcc0ytCHHY+qX8e2YTyYzj/55YwHrlYY6ZUs7ZW1Uzwk8SBIleCU5XSvphkrVjCZ38Dcw5ou+93lkhoaUZ+jBAdc5RbQku2rGCc7etJuMg7jkq/QxBL7QaS99qicTwdzuIyLa7Yi4gE4qQDEf1xUKKUk8C1LuAA4CH+6gsIiKECAg9+yDc8NvViQvmwO9+jH/cVwgffAKpAVikvtxS/GWf0XzaFLCkJYMBTamAbUdGqfUSBPmWDugl6Uy2dStaO7rwQZ4P/bTovR9kqG0zW1+Nb4NHJpOhpe3qDwpOpUj1JECNAXeZ2XNAu5n7zrkv9mqpRGTYijY3wj+vzp9559+I7nUYqdLK/i0U2UXq31iS5BtPL2LG8iSewVEblLH1yFH98vi+C7BQCDdl0+ySQWva4yCSGkcoIkNETwLU6bmbiAxRZkajC1GfMlKBoyJi1Prpfh2fZivqoNCi4ek01C2Gfg5QzYx36uHge+asSgsc3PlRI68sTvDU0ROo7qUtVAuJ1y/F/eSreN+4gOCu6+Htl7MZnoftdiDB579OcoC2PxUR6W3dDlCdcz83s/2BE4FRzrnDzGx7oKLPSici/cbzPGYnQpz+2AJeWJDd635MSYjf7TmKPceEidFPM7RD4c7zB2AMaoML8d1n8y/5PGtFilcWJTlgrNdngXyYAO7+O9TXEfz6POywE7EjT4Z0CkJhXKKVRDSuLUNFZMjodoBqZmcD3wSuAY7JJbcAVwK79H7RRKQ/LU6H2fP2WSxtXT22cEFzmuMfnMdjn53IDtXWL5MpgrJKvJFjYfH8jpmVNQSV1e2SwkE6Oyxg+VKIxZlSW5Xdick8Agce674feCIwXl1UeDH8B2c1ctD4qj4LUP1EK7z7evZOcyPu1r+2j0WnbIo/bWvSIe0EVEjSQtSnPdIOSn1HlZ/RzHWRItaTzfi+BeznnLuE1Yu2vAds3NuFEpH+5Xkej89pbhectvW9/yym0a37/vLd0VpSDudeCrE1dlKKROF7l9EaXz3OMp5qJXrn3+DsY+BHp8E5nyf53nTerAs489k6Tn5yKQ/MzbDcRbK76awlD6iKFn67nFQe7tPg3fkhqBlR+ICakTi/i5bnYcrMWBxEOfOZpWx840ym3TiT/e+ZxwtLs+OKRaQ49eQTpxz4NPf7ynfiMKzj6tciMuB83+fxT+sL5r+2qJWkM0p7ZQ/PzgVBQMvo9YhfcQu88hx88CZsMA123IuWskqCXCDoex7+fx+F+25ede6yL57LT91mXH3n3FVp933cyCbVER48cgLVln8b1K5UhzJ8Y6tqfvHS0g55BhwzpZxMpu/eCpPhKOHPfgmmv5b/gCNPzq5soBnZHSwLIux316fMWrF6E4P365IcePenPH3semxV0T89AyLSMz1pQX0GOG+NtG8AT/ZecURkIARBwEZVhVvgxpeFevRmsa4yQGNpJS17H0HrGT+mZf9jaCyrItNml/to0wq447rVJ5WW8eGUz3D1Bx2D0Pfqklz91vLsUkxrI8hw+maV7D2hpF2yZ3DDAWMZGe7bNUCdc6TX3xiOWmPBFDP4/FmkxkxUkJWH53n8b1Fru+B0JQd899lF/dYzICI905P/zLOB+8zsK0C5mb0PrAAO75OSiUi/yWQynLhRBRf9byn5lvP8/nY1VPvpfl/vMpMpHPhZJg31bfaF33In/jKn8FvaNe/Uc+bmlVR5axdMVpDg7/uNYk6z4+k5zdTEffYaX0KVnyZC3y9S3xqOEj38C4T3PQreew3Mh022IllSRnKwd+97IZoDiHrgB703Gc/zPO6fuaJg/v8WtpIIjNL+/PYlIt3Sk1n8881sB2AHsrtKfQq85JzTKHORIWBkOMNth4znpIfmkcisjlK/OK2CozYoJQiKazSPC4Wwqtrs5CggCEdoyBQeg9CcDtZ5knu5pZhWCptvGsc5RxCs3ZCBtZUIRUhU1ODttD/AoJ/kE5jPgpTP1W8t54UFraxfEeK729QwqQSivbBqhHOO9coLf8xVRz28fhi2IiI916O+jdw2py/lbiKyLjyfZWmf1owj5hs1oQwM4HaRYTLsOdrjnS9M5u2lSVYkA7YdFaMm7CgpwqHmiZIK4seeBtdcBoD33uucvEeaOz7Mf/wRG5RR5ge9shRTZy27PbVy8lZPuugHe2AK2dbNt1YY+975CcncF6KXFsC/Pmjgz/uO4ej1IoTXsWU6k8lw/EYVXPjS0rx/9m9sXUO1n0HNLCLFRx0bMmyYWXb5oSLQSISrpjezzS2z2PSmT9jqn7P49VtNNBAZ0HKFXMAIL8neo4wjJ/hMCCcpoeP4vWKQCQIyO+0DR38JQiFYNI/tkvPYqrZjHZaGjZ/sOIKwy7bKrcuM/t6SwWdpEOHZJXD/vAxzUxGae9ZmMKjVZXxOfWTequC0ra8/uYDlmd6pi1HhgOv2H9uhpXTvCSWctmkFbgC/FIpIYcPn3VCGrQw+dUF2m8rFLSm2HxVnTAzKbGACr5SFuPzV5fzu9dXjJxuSARe8tJT5TWku2rGqV7o318VgmXDTEo4RPvwLRPc7ClYsZ2wsxj0bj+H6D5r589vLaUoFHLZ+GefvWMu4cJomF2ZpEt5dlmRE3Gf98hA1fpr+bkJLm8/ziwOOe/BTWtKr6/rUaZVcuFMNZUXYYt3blifho/r8/4OpIDuxbdTIdZ9hHyHNoRPCvHvy+jzxaTPLEgH7TChhQolRTv8O0RCR7lOAKkNa2jxeWhpw9L8/pbVNS80B65Vwzb5jqBiAD6i6tMdVb9blzbtuej3nbFvDmGHyn+l5HisCn5QzfNZu8fSU55MqrVy1/eniGe/x7c034dRpZThnlPsBIZdkhYty1tOLuO/jxlXnVkc97jtiIpuX069B6pKUz1H/ns2ajYfXv1vPdqOinDIlOiS68TvTVdyZChzQOy3dETKM9jN8foPseriZzMB+ARSRrhVHf6dIH1mWDnHkfXPaBacAj8xu5o9v1eMGYO/ypa0B6QKxR+BgYXP3uhzNDN/3e7272vO8Ph8K4XkerRbmP0sch/17AVP/PpO975rLnbPTNLLuM9JdkKHKUlR7SUIuTeCF+N0by9sFpwB1iYAD757Nkkz/zYL3fZ9bZzR0CE5XuvjlZdT1Uvd2MauKGOPL8j9P32Dz2kivt+QHQdCr44dFpO8oQJUhy/M8Hp7dTKpAMHjVm3XUpfs/QI2HOg8oyyKd/1t6nseyIMJjCx2Xv9XEowscy4LoOgeqrYSYn45w88dJrpuRYFYiTHMvBIttZfBZHET571LHXZ+0cvA9c3h9cSvpINvde+qj8/n16ytI9HLnTl3a4+q38rdaN6UcLy1s6VZQ7nkeYTNC/tq/bsyM9+oKd+HPa0rnXeprqKkJpblm3zF5Z9H/YqeRVPkKJEWGs6H/NV2GLTNj5orCgUBDMiA9AIFATcSYVhPh3WUdy7Z+RZgR0cKBZiQSYWaLzz53zm63LWltzOfxoyeyfiy9Vq1OzYS44f1mfvDfxe3Sv7xZJb/YsYdjIj2fpWmfprQj7hvVoSDbionP/+ocR943kxsOHMcP/rMo7+m/e30ZZ2xRyehuxoCB+TQFHhFzxCx/UJMMsoFoITPqUtiEwsG453nEG+uxN16AV56F2tGw/9Ekq2pJhno2sS0IsmMg//F+/vU5txwRJWRDP0INgoDta3xePmEyv3x5KS8vbGVSeZgf7lDLFtU+kSKdnCci/UMBqgxZQRCw1/gSrng1f8vZlKoIUa//A4FKL8Vth4xn3ztnt+vOr4353H3Y+E4XxA/XjOXw++e2C04BlrZmOPr+uTz52QlUWs8n2HzSRIfgFOC6d+rZf2Iph473uzUmspEwN73XxMUvL6U+ERD24PObVPKzz9SSDuCI+2aSzDh8M5Yn8l8vcPB+XYoxXUyQCcxnftLnd6/X8dy8FsaWhvj+djVMGDmuw7ExD0aX+AWHT+w4Jlbw+ZkZJfVL4fzTV625CsBDtxH5yvdxO+9Pyg+tWiUiuz5q4boKgoC9J5RQG/M7/B0BLt1tJJV+pt83RRgIYTJsEMvwp91raM4YEc9RYhmcU3AqMtwpQJUhyznHFrURJlWE8251eNmuI6kegEDAOcd60TTPf2493lmW4s2lCTaribBlbYQRoXSnwU2jRQrOfP64PsXShKMy1rPyBF6IK19fVjD/0leWsfuYMZRZ5xUVeCGue6eRn7ywZFVaKoDrp9czsz7FxbuOXLWkUKiL3vTSsNHZHiCe5zG9wdjrjtVraE5fluDxT5v44XbVfKMqRKzNSgi1oTQ/32kEZzyxsMO11isPM606jHP5A/tIOgnXXd4+OF3pmsuIbrML4VAIb+4s+OBNGDUeNtmKlpLydluztjUylOLpY9fj1Efn8/LC1mxa3OeK3UexTbVP0Iu7KQ0GYZemMveaGCQLSIhIH1OAKkNajZfisc9O5OynFvLwrCYcMKYkxOW7j2TnkaEBCwSCIKDGS7LHSGOv0XGCIMC5ZJfBcksXYxKaC82+6kQqgLlNhcf7LWpOZ4dCdDHEdVnK49JX8gRxwNNzm0kGEPGNZMbxaUOaTaojecdilkc8JpeHoJNhBfUZn68+Pj/vGpoXv1LHyZtWMbbNu1sQBBw+KU7DbqP4xUtLaEhm62n3cSVcu98Yqr1kwcAo1NoMrz2fPzNeitfaApd/D+bPXp0eiRI///e0TJqaN0gNgoCJkRT3HDKG5SkjFTgqw0ZtKD3sglMRkXwUoMqQ5pxjpJfg+n1Gsjw9imTgKA8VTyDgnOvRrOLaqOEbeWeA+waj4j70cPedmOc4YL0Snp3bnDd/9/FxSn3X5Q5MK1Ku03GejckMf9t/LBHfozLicfW+Yzn83k9XBYuQbVm95aBxubVJO3mstPH20vxLhDng+fktHDsp3K41upQUp28c4+gNJ7E8GRD3ParDjhJLdj5uN5Mp2Kxn+x1FcPMf2genAMkE/PKbxH7zL5pKKvKX0zlKSVHaZujrcOjWFxHpDgWoMizESLVbW3SwBgJeUx1nbVnNlW90HFd75pZVVIV6vpWnCzKcsFEFv3p1GfVrjAsNe/CjHUYQLtD93VZJyDDyP/zpm1fRmHJ877lFzG3MfjE4ddMK/nPcZO75qIHXFreyRU2EEzeuYGQ4g7l1m8FdaEEDCzLUehlq2wyD6KpLOYiV4E/cAD79uOP1NtsWd9/N+U9sbcFmf4RN23bQbHwgIlIstMyUyCDSVLeE721bxSW7jqQmlp3mXh31uHiXkZy3bTURt3atwqP8JM8eO4k9xpesSttyRJSnjpnE+Gj3gsWqsOPASaUd0idXhNlvYiknPDh3VXAKcP30FbwzcwHfL53Dvxoe4Mcv/p4p/7yMyuWL8LuIsivDsMWIaN48A3YeE++1he4TsRL42g8h3zJU4WjnEW79sqLYVlVEZLBRC6rIIFNGkjM3iXPclPVIBBD1oDaUwQVrvz2mc45J0RS3HDCKFWkjcI6KMJRb95etirkUV+01msPvm8v0Zau738/asprL8oxNvWDLGAf/7xa8R/7VPuPpB4hfeA1NYycVfOxyS/HXfcew5+2zSawx3uHHO1RnW5J7SRAEtI5fn9hlN8JNv4d3X4fKGjjiZNyYCVBRDSvyrxTBBpsU3Y5QgRciPnI8rYSJ9+DvKyLSnxSgigxCLshQ42VW9YH0xi6dzjlKSFESapvWs2uM8JM8eMQ4ZjdmeHNJgskVYTaqjvK959qveRry4KTRKcqu/FfHi6SS8OdfEv/R72gJx/IGUM45Ni4NeP2kyfz+jTqemdvChLIQ525Xw1jX0G4Gf29Iez5NoyYQ+cYF+KkEzjyS8TLAEf/C1+GPF3Qs4+bbEVTV9mo51kVgHgtSIX77Wh1PzW1mdEmC87arYevaMKVac1REiowCVBHpNc45Ki3JFuWwVWV2q8q6IGBMaYj5TauDxmk1UUrfe7HwhT5+j9CcTyjFkRk3iZZIvMMh5gLGhpJctEMlTdtWEvEgSoYPPpgHVVP75Lkl/DD4uVlNuZbR1La7Ef7mBdnW1aWLIBKD/Y6Eo06hJdzDNb/6iJnxYZPH7rd/sqrF+YO6JM/ObeY729Rw7lblxG3gJw2KiKykALUbzIwkPsnAKPUDCLQFn0hXVnZt1/hpzt22hu88u7oVNXDkH9PZVksjXPJd/B32JP6V82iJdgxSASxIU5abnTUQndWJcJT09nsS3XRbLJmAUIhkSTmprtbl6kcNLsTXnljQYTgEwBWvLeO0zSoZ37u72oqIrBNNkupCE2H+V2ec/PgSDr1/ARe+1sjCTBRMVSfSHUEQ8LkppRw/tXxV2nvLEjRuvF3hk6Zujpv1Yfb3/z2N//G7RT3ZKBMENMdKaaqooamkoqiCU4CGlPHa4taC+U/Nacbr6guDiEg/UgtqJ1oIceWbK7j0ldW77Ly+uJU/vVXHc5+bxAYxpwkGIt1QTpLf7FbLD3eo5eVFrVRGPCprHBz3Vbj1L+0PjsXxjvsKwZ8uXJ12/z+JbLwlCU9vWWuli3jZL6J42syIJVrwWpsAyMRKSUTjeq8VGWb0bt+JpUmvXXC6UlPKcdaTC7n9oNGUaHKBFJGVrYzF+GFeQopJUVh/Uhjnsl/uUgccQ3jLHeCeG2HZEmyjzbHP7EVw81WwbPHqkxOt2TGfauRbKxUhx/ajY6u2VV3THhNKCNZhFYje4uOIL/gU/nQhzHwfgNCG0wid8SNaRk3o4RYUIjKY6e2+AM/zeGJO/p11AP47v4WGdBE1O8iw1kqIeakI983N8NhCx+IgSsr8gS5WXtltXbMBdCIcpWm9jUh//WfY507DLVlA8IuzYcY77U/a7QDSkeKYcDQYlVmaP+8zhpJQx/esH+1QS3UvLsu1LuIr6uBHp60KTgH46F340enEGwos5SUiQ5JaUDsRdNEKVXxtVDIcNRLmkleX84c3lq9KC3tw7X5jOWh8mEiRtzs552jFIzRpKrz/VsdJiCPG4Hbcm0yRrSc6mDjn2DCe4bWT1udPb9Xx1JwWxpb6nLNtDdMqfeJF0BMUxsF9N2eXGVtTshUeupXQ577WywuIiUixUgtqAUEQsM/EkoL5O4yOUR5SiCoDy/M8np6XaBecAqQCOOWR+cxPDJ5/8eayKrjketj/s1BSBmUVcNiJcOFfaS4p7+r0ohMyKGlpoLRhGSWtTfgDPQnJBYz2E/x02wpu3jXK3/cZwfZVrteGKfm+TygUWuvJbH6iFaa/WviAd17BTxWe6CUiQ4taUDsxIuL4+lZVHT78Y77xx71HU2bpHi9kLtKb6jM+F/1vUd48B/zl7Xou3rGCIFP87U5BENBYWkn4C2cTPvY0wEjGSrMtZoPsH60k2YL30G3wwL+gpQmrHkn8xDNIb7sbreH8W7T2FwvSNC+ZT7i6rFd6gZKEWJLyuOv9RuY3pTl4Uimb10So8lI9GgvtQiGoHgFzZuY/oGokztdaWCLDhQLUTpSQ5rxtqzhwvVIufaWORS1p9h5fwre2qWZMWFsEysBLO2NuU+EWsA/rk6QcFOdo1PxSeKRiZQNdjG4zMyLpJF4qSRAKY4B33a/g+cdXH1S3GP54AaFTvk1o3yNJD5G3jiQ+989J8uVH568Kdv/wRh0bV0d46MgJVFui0/PbXSsUIfzZU+Ct/+U/4LNfJOWHBt2XFRFZOwpQu1BGij1GGjscNIpkAGV+gBckNQBVikLMC9huVIwnPs0/oW/vCSWECdDozb4RDtJEF82DW/8Msz7EG7se/qnfJmgbnLZ161+I7rQP6UE4ZCGfpWmfLz86u8Pb4ft1SS54aSmX7VRFyHWv9d45R3q9qYSO+iLc/ff2mceeTnrcZDUKiAwjClC7wTlHjBQxA33SSzEpsQwX7jySXT+d1SFIqIh4HDOlnCDofiuWdJ9vRnT6a3D5uavSLJPBvfdG4ZNamrCmBhgCAarv+zzwYVPB7+o3vb+CH+1QQ20Pht62hqNED/8C4X2OyI5HNQ823YZkSQXJkLr3RYaTwTODQkQ6cM4xpRTuPXwCE8pWf9/cemSMp45Zj1GhgZ+dPVRFmxvg6ovaJ7a2ZCd3dSYS6btC9bPFLYVXiEhmHHl2Vu1SIhShsXIEzbseRPMuB9BYUavgVGQYUguqyCAXIc2eozyeO3Yi9UmH70FVCMq9NC5Ql2hf8RrroWF5+8SWJvD9bJDauKLDOW6TrcgMgdZTyE5qO2C9Ei55eWne/M1ro8S9tX/9BVpWTGRYUwuqyBAQBAFVlmRSNMWEcIoy69kMalkLBZZTcndej3fWTyAWb59ROwr7+s9IDJENB5xzTKkIsd2ojs/HgN/tOYpKv7jX4BWR4qUWVBGRtRCUVeLVjoKlayzz9fF7BP/+B/brf+JmvA3zZsHUzQnW25CWePmQ+uJQYUnuPHQcv3m9jr++vZymlGObkTF+s8coNq2AYM1NF0REukkBqois4nkedRmfuiSkAkdN1KMmlMEUaHTQGi+j5JsXwC/OgnSbmepmcNBxtJRWEmy3B7a9rd7edQgFpytVkl38/5tbVZEOIO47yk3L8InIulGAKiIAmOfxQZPPiQ/NY8by7HaT5RGPS3YZyWfXj1GiTSbbCZyjdb2NiP3mX/Dw7fDhdBg/GQ49kUT1SNJmMEzGUVqQptpYteCuYlMRWVcKUEUEgMXpMHvf+QkNydVBVUMy4KynFjKhfAJ7j/I0cWUNaTMaK2oJHX8GfjJBEI6QUnAmfcDM1Cotw4omSYlIdk3LTxrbBadt/frVZaTTAdFUgtBa7rU+lKUDRyKk4FR6XzzZQtmSeZS+/RJl82ZS0tqE6X9QhgG1oIoIZsZz81rz5v1++zhHhuZT8YcroH4p4a13hn2OoLmsikAtOiJ9prS1CX7zQ2iz+YM3YgylP/49zTWj1aMhQ5paUEUE5xxbjYh2SP/J5lG+8N7djLvk/+CVZ7PjLG+/Fs75PCXLFqolR6SPRDMp7LpfY2vuTLZkAVzwdWLNHdfZFRlKFKCKCJlMhmOmlBFu844Q8uCUMSnK//33jic0N+KuuYxIStuoivSFUHMjvPRk/swlC/HWXN5MZIhRgCoiAIwKZ7j/yIlURbNvC1uOiFE+/cWCx9vbLxNqbe6v4okMK5Zo6Xw5hKWL1IMhQ5rGoIoIAJ7LsGO1x6snTGJOU4aQQdVLXZykMagifcLFS7FwBFLJ/AeMnahZ/TKkqQVVRFZzATVeki3LM2xalsHbaqfCh260BZlYST8WTmT4SJZUwMHH5c1zUzYjU1nTzyUS6V/9EqCa2UQze9LM3jWzd8zsm7n0GjN71Mxm5H5WtznnB2b2oZm9b2YH9kc5RYYaz/Ood2E+SYT5qDXEsiCCed3/t09VVMP+n+2YEYliXz1vyOwrL1JsUmYEh38BDjsJQuFsohlstxuceyktkfjAFlCkj/VXF38a+K5z7lUzKwdeMbNHgVOBx51zl5jZecB5wPfNbFPgBGAzYBzwmJlt5JzTfosi3eTM450Gjy88PJeP6lMA1MZ8/rDXaPYZGybajZ2hEuEo3vFfw992V7jrBlixHLbYHo44mZaKmgHpYlw57k7dmzLUNUdihI89ncghx2PNTRCLkyopJeFHBrpoIn2uXwJU59x8YH7u9wYzexcYDxwJ7JU77AbgKeD7ufRbnHMJYKaZfQjsCDzfH+UVGWzSFiIRQNwHL8gGnovSIfa58xNa0qsDuaWtGU58aB7PHLseW1V0b2ealnAMb/MdCW+0JV4mTToaJ0X/T85oIURdyuPlRa3EfGPrkVGq/AwR9L1Vhq6UeaRKK6G0cqCLItKv+n2SlJlNBrYBXgRG54JXnHPzzWxU7rDxwAttTpuTSxORNloJMbMJLn91KR8uT7HdqBjf3qaGCXHHP6Y3tAtO2zr/+SXccsAoSkh163GCICARikBoYFpuGolw2WvLufL1ulVpIQ+u3nsMh68X7VZrsIiIDB7Wn91kZlYGPA1c5Jy708yWO+eq2uTXOeeqzewq4Hnn3E259GuBB5xzd6w8tr6+flXBZ8yY0W/PQaRYxCuqeK6lkjOeXNwu3Tf49xHjuefjJq5+a3nec6uiHs8dPoKWBZ/0fUHXUSwW4y1vLCc9sjBv/ovHTiC06EPtqiMiMohMnTp11e+VlZUduuX6rQXVzMLAHcDNzrk7c8kLzWxsrvV0LLBy5eE5wMQ2p08A5hW6dtsnWaxmzJgxKMo5kFRHXWtbR0uCCN+86ZMOx2QcnP74Qm44YFzBAHVyRYSq8hImlBd/fTcS5tJ/zy+Y/9d3G/n1zhvjMtlWVL2OuqY66prqqHtUT11THa2d/prFb8C1wLvOuSvaZN0LnJL7/RTgnjbpJ5hZ1MzWB6YCXa3IKDKsfFyfJpHJ3wMytzFNRcTHKzBU9Cc71lLhDY6xm6nAmNtYuAt/5ooUSTWeiogMKf21DuquwMnAPmb2eu52CHAJsL+ZzQD2z93HOfcOcCswHXgIOEsz+EXac3Q+PMc3x22HjKckZG3SssHpTqPCg6ZLvMQP2GFU4eWs9poQJ+ppRr+IyFDSX7P4n4OC0373LXDORcBFfVYokUFuw8oIEd9I5mlFHVsaojZqrD/a563PT2bmijStGcdGVRGq/DSRbk6OKgYxMvx8pxE8NKupQ0heHvE4YWoFQSYxIGUTEZG+oZ2kRAapKj/Nb3Yf1SHdM7hmvzHU+Gk8l6HWS7J9VcDuI2CUnxh0yzI555hc4njgyAlMqgivSt96ZIynjlmPkaHBE2yLiEj39PsyUyLSO0Iuw9Hrx9hixHpc+vIyPqxPsu3IKN/bvpYJ0YAgaB+IDuaF7SNk2HWExzNHT6A+6fA9qAxDuaUH9fMSEZH8FKCKDGJx0mxVYfxtn1paA6PEd4SCodmiGAQBlZakMro6TbGpiMjQpABVZJBzzhEhTcSAwTHvSUREpFMKUEV6wDyf5WmPAKPcDwg57WAkIiLS2xSginRTPVH+8W4DV7+9nMZkwOHrl3He9jWMCafBqelSRESktyhAFemGehfhmAfn8fLC1lVp179bzx0fNfDi8ZMYF0oOYOkGn4wXojljRDxHjIwmOomISDtaZkqkC2bGG0uT7YLTlRqSAb94cQkp03e97sjgMzcV4fsvLOeAe+dx4iOLeXGZ0Uy465NFRGTYUIAq0gXP87nxvRUF8+/5uJHGjP6VuuJ5Hu81Gtv88xOueaeed5cleWJOMwfc/SlXv9NAqzp0REQkR5+qIl0wg9Jw4X+VmG9YoX3SZJX6jM9XHl+Qd+erC15aSl1ab0ciIpKlTwSRLqTTaU7frLJg/pc3q6LS1ySprtSnYPqy/FuSOuDFBa14nt6SREREAar0spDnEXEBoSEWaEwu8/hyniB1SlWEr29ZhQVabqorXbUyqxFaRERW0qAv6RUhFxBrqIPH74FPPoD1N4F9jqC1vJr0EOj/LiPFL3as5uRNKrnqzToakgGf37iC3cfFqfGS2tGoGypCsHltlLeXdmxFNeAzY2IEgVZDEBERBajSC3yD2Mz34YKzIZNrSXz9BbjvJmI//gPN608jGAIRXBkptq00rtmjhgwQJiAIEgpOu6nCS3PNfmPY8/bZJNYYh/rznUZQHdIwCRERyRpa/bAyIKLNjXDFD1cHpyul03DFD4m1NAxMwfqAcw6CNH6QJggUUPWEc46NSwNeP2ky/7dlFVuMiHLI5FKePHo9Tp9WRhQNkxARkSy1oMo681bUwYq6/Jn1y/BWLIdYWb+WSYqTuYCxoSQX7VBJU6aSqAcRMjiXGeiiiYhIEVGAKusu6CK4yAy/4MPMqAvCLG11pJ1jRMyjJpTBuqqrYcKCNGUGuOwMfhERkbYUoMo6c5U1WCQGyY47LRGN4Sqr+79QAygwn7dXwIkPfcrcxmy3dVnY44rdR3H4pCjxfuzKDoVCOMAFgYYkiIjIoKExqEOE53lEUwliieZ+/9aRiJfBF7+RP/OUb5GIl/dvgQbYorTP/netDk4BGlMBX31iAW8tD7BurGpglt3+c3EmwrIgQrqHW6m2EmJWIsyPX17BqU8u465P0yxz0W49toiIyEBTC+oQEE+24H/wJtxzEzTWE9p6Zzj0BFrKq8mzaU+vS5tHaqd9CY+bBP/8E8ybBeMmwUn/R2q9KUNimanu8n2fm99Z0WGW+ko/fn4Jdx08mhJSBa+RxuODRo9vP7OQFxe0EPGNkzeu4Ic71hIKdf0vm8TnjpmtnPXUwlVpd3zYwPiyEE8dvR4j/fyL5a8UMiPatAKSCYhESJRUaPqSiIj0KwWog1ws1Yr/99/Bsw+tTpw3G564j/ilN9BUPSo787yPJcJRklM2J/L9X+OnU2RCYZKRWL88djFJOuN/iwoHgO8uS5AIjJICfRdmxsdNHrvfNmvVl4tkxnHt9HqemtvMvw+a1GUZlqZ9vt4mOF1pbmOa7/9nEVfvWUvE5Q8546lW/MfuhntvguZGiMWJHXwcwSEn0ByJd/nYIiIivUFd/INcaNni9sHpSq3NcP0VRFL9t/C5c45EOEpzvIxEODrsglOAsDm2rI0UzJ9SFSFsheul2YX43n8W5235/qg+xRvLg063A/U8j8dmNxeceHT3x43UF9jzPhxk8O++AW65OhucArS2wF034N30B6Lpwq2+IiIivUkB6iDm+z68+GThA15/gVCypf8KJASZDKdOqyRU4D/rZ58ZQYVfeCZ/a2A8M7e5YP69nzR3GqCaGcuTha8fOEgXiF6jzY3w4G35M595gHBLY8HrioiI9CYFqINdZ3veD6Oxn8VkdCTNfYdPoDq6+m8T8Y1Ldx3J9iNCnc6mN4PKSOG/6ei43+ljZzIZ9ptYWjB/89ooZX6BCLWxvuNmCys5B8sWd/rYIiIivUVjUAexTCYDO+0Dt12T/4DtdicdK+nfQgm+C9i51uOVEyaxsCUgkXGML/WpCmUIu867yav8DGdtVc2FLy3Nm3/S1FLS6c6nLE0oMQ5Yr4RHZrdvifUMfr/XaCr9DHlj5Gis0+tSqs0WRESkf6gFdZBLVdbCAcd0zCirgC9+g4Sn7yADwgXUeEmmlabZuiLDSD9JuDu7JQUZvrJZJTuN6Rgs/mr3UcRb67u8RLmluGbfMVyy60jGlIQIebDXhBL+e9wkNq+gYAtuprQcJk3Nf9FR4wjKKrsuv4iISC9Q9DLIJcJRvM99BX+nfbIzrxvqYfvdYK/DaC6tzHbNyqBSSYJbDxrLhyvS3Dezkeqoz9FTyhkRDvh0xicwoutAsYIEZ24S54QpZTiDqDlKLd3pxLXWaAml51wKPzsDli5qU6Bq+OFvaS0pJ3/Tq4iISO9SgDoEtERi2JTNCX/zQiyTJhONkQ6cgtNBrMKSbFdl7Lh9Bc45Mpmer8bgggxV3upW265eDs45miprif/yOry5s2D2hzB+Mm7ihrSUVmgnKhER6TcKUIcI5xxJzwfPz07VlkHPOdfleNO+eMzmWBk2ZXNs6hY457KtrgpORUSkHylAFZEOVgWmIiIiA0CTpERERESkqChAFREREZGiogBVRERERIqKAlQRERERKSoKUEVERESkqChAFREREZGiogBVRERERIqKAlQRERERKSoKUEVERESkqChAFREREZGiogBVRERERIqKAlQRERERKSoKUEVERESkqIQGugAiA8F5HnXpEEsSAQAjYx7VXhpcMMAlExEREQWoMuwkCfH0/DRffXwWdbkAtSbmc+1+Y9htlE/YZQa4hCIiIsObuvhl2JnVAsc9MHdVcAqwrDXD0f+ey6ct+pcQEREZaPo0lmElbSEufXkZLk+eA3796jLSpo4FERGRgaQAVYaV5sB4e2miYP5bSxO0ahiqiIjIgFKAKsNK3HNsUhMpmL9JTYSY348FEhERkQ4UoMqwEnZpztu+tmD+udvWEArS/VgiERERWZMCVBl21i+FGw8YS1l49cu/POLxz4PGMSmeb3SqiIiI9CfNBpFhJ+rSHDohzBsnTWJBc4ABo0s8qv0MnpaYEhERGXAKUGVY8lyGWi9DbVmbRDWeioiIFAV18YuIiIhIUVGAKiIiIiJFRQGqyCDheR4VNbV4vkbmiIjI0KZPOpEiZ2bUBWGen5vgHx/4VEXq+OrmVWxQ5lFqqYEunoiISK9TgCpS5JYFEQ66Zw7v1yVXpd343grO3KKK87evohQFqcOZD8RaGrDlSyGTwdWMJBEvI23qIBORwauoA1QzOwj4Hdn34Gucc5cMcJFE+pXzfK56bXm74HSlP721nJOnVbJp6QAUTIpCyGWIzXgHfvtDaGoEwMIRYl/8BqldDiARKrxrmohIMSvar9hm5gNXAQcDmwInmtmmA1sqkf5Vn/G59p3lBfOvfWc5oVBRf8+UPhSrXwYXf3tVcApAKgnX/orw7A8xs4ErnIjIOijaABXYEfjQOfexcy4J3AIcOcBlEulXDmhOF16gdUUywGn91mEpZAaP3AFBgc0lbvkzkWRr/xZKRKSXFHPTy3jg0zb35wCfyXfgjBkz+qVA62qwlHMgqY7aK6sZyeGTS7n9o8a8+SdOKeXjjz8imew4BGA4Gw6vozGVFYyd/VHhAxbMpqV+OR/VrcibPRzqaF2pjrpH9dQ11VFHU6dO7TS/mAPUfH1TeduKunqSxWDGjBmDopwDSXWU3893jvDArKYOLanbjIyy7eg4VTZpgEpWnIbL6yjkGWwwDd76X958N2ED4lXVTB0xukPecKmjdaE66h7VU9dUR2unmLv45wAT29yfAMwboLKIDJjxkTQvnzCZEzeqoDziMaYkxM93GsFdh42nytRyOlylAwf7HQUFxiDbiWeQDEf7t1AiIr2kmFtQ/wdMNbP1gbnACcBJA1skkQHgAsaFk1y5WzXnbx2ntKSEaj+DCxIDXTIZYK3l1cR+chX85kdQtySbGC+F084lOW4yTgOURWSQKtoA1TmXNrOvAw+TXWbqOufcOwNcLJEBE3ZpmufNZPzUqbhgoEsjxSBtRvP604hdfD1eYz0EAa6iikRJOZ3MrRMRKXpFG6ACOOceAB4Y6HKIiBSrIAhojpdBvGx1ooJTERnkinkMqoiIiIgMQwpQRURERKSoKEAVERERkaKiAFVEREREiooCVBEREREpKgpQRURERKSoKEAVERERkaKiAFVEREREiooN1q3w6uvrB2fBRURERGSVyspKWzNNLagiIiIiUlQUoIqIiIhIURm0XfwiIiIiMjSpBVVEREREiooC1F5gZhPN7Ekze9fM3jGzb+bSa8zsUTObkftZ3eacH5jZh2b2vpkdOHCl719m5pvZa2b279x91VEbZlZlZreb2Xu519POqqP2zOzbuf+zt83sn2YWUx2BmV1nZovM7O02aT2uFzPbzszeyuVdaWYdJi8MVgXq6PLc/9ubZnaXmVW1yVMdtc87x8ycmY1ok6Y6Wp1+dq4e3jGzy9qkD7s66hXOOd3W8QaMBbbN/V4OfABsClwGnJdLPw+4NPf7psAbQBRYH/gI8Af6efRTXX0H+Afw79x91VH7+rkBOD33ewSoUh21q5/xwEwgnrt/K3Cq6sgB7AFsC7zdJq3H9QK8BOwMGPAgcPBAP7c+rqMDgFDu90tVRx3rKJc+EXgYmAWMUB11eB3tDTwGRHP3Rw3nOuqNm1pQe4Fzbr5z7tXc7w3Au2Q/SI8kG3CQ+3lU7vcjgVuccwnn3EzgQ2DHfi30ADCzCcChwDVtklVHOWZWQfaN71oA51zSObcc1dGaQkDczEJACTAP1RHOuWeAZWsk96hezGwsUOGce95lP0H/3uacQS9fHTnnHnHOpXN3XwAm5H5XHbX3G+B7QNuJK6qj1c4ELnHOJXLHLMqlD8s66g0KUHuZmU0GtgFeBEY75+ZDNogFRuUOGw982ua0Obm0oe63ZN/ggjZpqqPVNgAWA3/LDYO4xsxKUR2t4pybC/wKmA3MB+qdc4+gOiqkp/UyPvf7munDxZfJtmSB6mgVMzsCmOuce2ONLNXRahsBu5vZi2b2tJntkEtXHa0lBai9yMzKgDuAbznnVnR2aJ60Ib2cgpkdBixyzr3S3VPypA3pOiLbMrgt8Cfn3DZAE9lu2UKGXR3lxlAeSbarbBxQamZf6OyUPGlDuo66qVC9DNv6MrMfAWng5pVJeQ4bdnVkZiXAj4Cf5MvOkzbs6ignBFQDOwHnArfmxpSqjtaSAtReYmZhssHpzc65O3PJC3PN+OR+rmzyn0N2PM9KE8h2Uw5luwJHmNknwC3APmZ2E6qjtuYAc5xzL+bu3042YFUdrbYfMNM5t9g5lwLuBHZBdVRIT+tlDqu7uNumD2lmdgpwGPD5XHcrqI5W2pDsF8I3cu/fE4BXzWwMqqO25gB3uqyXyPYUjkB1tNYUoPaC3Leka4F3nXNXtMm6Fzgl9/spwD1t0k8ws6iZrQ9MJTtYeshyzv3AOTfBOTcZOAF4wjn3BVRHqzjnFgCfmtnGuaR9gemojtqaDexkZiW5/7t9yY75Vh3l16N6yQ0DaDCznXL1+8U25wxJZnYQ8H3gCOdcc5ss1RHgnHvLOTfKOTc59/49h+yk4AWojtq6G9gHwMw2IjvJdQmqo7U30LO0hsIN2I1s0/ybwOu52yFALfA4MCP3s6bNOT8iO5vvfYbZzD1gL1bP4lcdta+brYGXc6+lu8l2GamO2tfRz4H3gLeBG8nOjh32dQT8k+y43BTZIOK0takXYPtc3X4E/IHchi5D4Vagjj4kO0Zw5Xv31aqj9nW0Rv4n5Gbxq47avY4iwE255/wqsM9wrqPeuGknKREREREpKuriFxEREZGiogBVRERERIqKAlQRERERKSoKUEVERESkqChAFREREZGiogBVRGSAmdmDucXiC+VfbWY/7s8yiYgMJC0zJSLSS8yssc3dEiABZHL3v+acu7njWR2ucSpwunNut94voYjI4BAa6AKIiAwVzrmylb/ntoU83Tn32MCVSERkcFIXv4hIHzOzn5nZTW3uTzYzZ2ah3P2nzOx0M5sGXA3sbGaNZrY8l3+9mV3Y5vyvmNmHZrbMzO41s3Ft8pyZnWFmM8yszsyuym2lKCIyaChAFREpEs65d4EzgOedc2XOuao1jzGzfYCLgeOAscAs4JY1DjsM2AHYKnfcgX1YbBGRXqcAVURkcPk8cJ1z7lXnXAL4AdkW18ltjrnEObfcOTcbeBLYuv+LKSKy9hSgiogMLuPItpoC4JxrBJYC49scs6DN781AGSIig4gCVBGRvtdEdlb/SmM6ObarpVXmAZNW3jGzUqAWmLvWpRMRKTIKUEVE+t7rwB5mtp6ZVZLtli9kITDBzCIF8v8BfMnMtjazKPBL4EXn3Ce9WWARkYGkAFVEpI855x4F/gW8CbwC/LuTw58A3gEWmNmSPNd6HPgxcAcwH9gQOKG3yywiMpC0UL+IiIiIFBW1oIqIiIhIUVGAKiIiIiJFRQGqiIiIiBQVBagiIiIiUlQUoIqIiIhIUVGAKiIiIiJFRQGqiIiIiBQVBagiIiIiUlQUoIqIiIhIUfl/fS+bnYr/1YEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "sns.scatterplot(x=\"Tuition\", y=\"enem_score\", hue=\"Tablet\", data=data, s=70).set_title('ENEM score by Tuition Cost')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we understand the problem let's look at the solution. We can also say what would be necessary to make association equal to causation. **If $E[Y_0|T=0] = E[Y_0|T=1]$, then, association IS CAUSATION!** Understanding this is not just remembering the equation. There is a strong intuitive argument here. To say that $E[Y_0|T=0] = E[Y_0|T=1]$ is to say that treatment and control group are comparable before the treatment. Or, when the treated had not been treated, if we could observe its $Y_0$, its outcome would be the same as the untreated. Mathematically, the bias term would vanish:\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = E[Y_1 - Y_0|T=1] = ATT\n",
    "$\n",
    "\n",
    "Also, if the treated and the untreated only differ on the treatment itself, then, $E[Y_0|T=0] = E[Y_0|T=1]$ and we have that the causal impact on the treated is the same as in the untreated (because they are very similar).\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "E[Y_1 - Y_0|T=1] &= E[Y_1|T=1] - E[Y_0|T=1] \\\\\n",
    "&= E[Y_1|T=1] - E[Y_0|T=0] \\\\\n",
    "&= E[Y|T=1] - E[Y|T=0]\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "In this case, the **difference in means BECOMES the causal effect**:\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = ATT\n",
    "$\n",
    "\n",
    "Additionally, if the treated and the untreated only differ on the treatment itself, we also have $E[Y_1|T=0] = E[Y_1|T=1]$, that is, we make sure that both treated and control groups respond similarly to the treatment. Now, besides being exchangeable prior to the treatment, treated and untreated are also exchangeable **after the treatment**. In this case, $E[Y_1 - Y_0|T=1]=E[Y_1 - Y_0|T=0]$ and\n",
    " \n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = ATT = ATE\n",
    "$\n",
    "\n",
    "Once again, this is so important that I think it is worth going over it again, now with pretty pictures. If we make a simple average comparison between the treatment and the untreated group, this is what we get (blue dots didn't receive the treatment, that is, the tablet):\n",
    "\n",
    "![img](./data/img/intro/anatomy1.png)\n",
    "\n",
    "Notice how the difference in outcomes between the two groups can have two causes:\n",
    "\n",
    "1. The treatment effect. The increase in test scores is caused by giving kids tablets.\n",
    "2. Some of the differences in test scores can be due to tuition prices for better education. In this case, treated and untreated differ because the treated have a much higher tuition price. Other differences between the treatment and untreated are NOT the treatment itself. \n",
    "\n",
    "The individual treatment effect is the difference between the unit's outcome and another theoretical outcome that the same unit would have if it got the alternative treatment. The actual treatment effect can only be obtained if we have godlike powers to observe the potential outcome, like the left figure below. These are the counterfactual outcomes and are denoted in a light color.\n",
    "\n",
    "![img](./data/img/intro/anatomy2.png)\n",
    "\n",
    "In the right plot, we depicted the bias that we’ve talked about before. We get the bias if we set everyone to not receive the treatment. In this case, we are only left with the $T_0$ potential outcome. Then, we see how the treated and untreated groups differ. If they do, something other than the treatment is causing the treated and untreated to be different. This something is the bias and is what shadows the actual treatment effect.\n",
    "\n",
    "Now, contrast this with a hypothetical situation where there is no bias. Suppose that tablets are randomly assigned to schools. In this situation, rich and poor schools have the same chance of receiving the treatment. Treatment would be well distributed across the tuition spectrum.\n",
    "\n",
    "![img](./data/img/intro/anatomy3.png)\n",
    "\n",
    "In this case, the difference in the outcome between treated and untreated IS the average causal effect. This happens because there is no other source of difference between treatment and untreated other than the treatment itself. All the differences we see must be attributed to it. Another way to say this is that there is no bias.\n",
    "\n",
    "![img](./data/img/intro/anatomy4.png)\n",
    "\n",
    "If we set everyone to not receive the treatment so that we only observe the $Y_0$s, we would find no difference between the treated and untreated groups. \n",
    "\n",
    "This is the herculean task causal inference is all about. It is about finding clever ways of removing bias and making the treated and the untreated comparable so that all the difference we see is only the average treatment effect. Ultimately, causal inference is about figuring out how the world works, stripped of all delusions and misinterpretations. And now that we understand this, we can move forward to mastering some of the most powerful methods to remove bias, the weapons of the Brave and True, to identify the causal effect. \n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "So far, we've seen that association is not causation. Most importantly, we've seen precisely why it isn't and how we can make association be causation. We've also introduced the potential outcome notation as a way to wrap our heads around causal reasoning. With it, we saw statistics as two possible realities: one in which the treatment is given and another in which it is not. But unfortunately, we can only measure one of them, where the fundamental problem of causal inference lies.\n",
    "\n",
    "We will see some basic techniques to estimate the causal effect, starting with the golden standard of a randomized trial. I'll also review some statistical concepts as we go. I'll end with a quote often used in causal inference classes, taken from a kung-fu series:\n",
    "\n",
    "> 'What happens in a man's life is already written. A man must move through life as his destiny wills.' -Caine  \n",
    "'Yes, yet each man is free to live as he chooses. Though they seem opposite, both are true.' -Old Man\n",
    "\n",
    "## References\n",
    "\n",
    "I like to think of this book as a tribute to Joshua Angrist, Alberto Abadie and Christopher Walters for their amazing Econometrics class. Most of the ideas here are taken from their classes at the American Economic Association. Watching them is what's keeping me sane during this tough year of 2020.\n",
    "* [Cross-Section Econometrics](https://www.aeaweb.org/conference/cont-ed/2017-webcasts)\n",
    "* [Mastering Mostly Harmless Econometrics](https://www.aeaweb.org/conference/cont-ed/2020-webcasts)\n",
    "\n",
    "I'll also like to reference the amazing books from Angrist. They have shown me that Econometrics, or 'Metrics as they call it, is not only extremely useful but also profoundly fun.\n",
    "\n",
    "* [Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com/)\n",
    "* [Mastering 'Metrics](https://www.masteringmetrics.com/)\n",
    "\n",
    "My final reference is Miguel Hernan and Jamie Robins' book. It has been my trustworthy companion in the most thorny causal questions I had to answer.\n",
    "\n",
    "* [Causal Inference Book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)\n",
    "\n",
    "The beer analogy was taken from the awesome [Stock Series](https://jlcollinsnh.com/2012/04/15/stocks-part-1-theres-a-major-market-crash-coming-and-dr-lo-cant-save-you/), by JL Colins. This is an absolute must read for all of those wanting to learn how to productively invest their money.\n",
    "\n",
    "![img](./data/img/poetry.png)\n",
    "\n",
    "## Contribute\n",
    "\n",
    "Causal Inference for the Brave and True is an open-source material on causal inference, the statistics of science. Its goal is to be accessible monetarily and intellectually. It uses only free software based on Python.\n",
    "If you found this book valuable and want to support it, please go to [Patreon](https://www.patreon.com/causal_inference_for_the_brave_and_true). If you are not ready to contribute financially, you can also help by fixing typos, suggesting edits, or giving feedback on passages you didn't understand. Go to the book's repository and [open an issue](https://github.com/matheusfacure/python-causality-handbook/issues). Finally, if you liked this content, please share it with others who might find it helpful and give it a [star on GitHub](https://github.com/matheusfacure/python-causality-handbook/stargazers)."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
